Solution 1:

Your application of the mean value theorem is not always valid as there is no reason why all the points on the line between $(a_1,a_2)$ and $(a_1,b_2)$ has to be in $D$. To close this loophole you need to use the missing ingredient. If $a,b\in D$ then since $D$ is connected we know there exists a smooth curve (see e.g. this answer) $\gamma(t):[0,1]\to D$ such that $\gamma(0) = a$ and $\gamma(1) = b$. You can now apply the mean value theorem to the function $g:[0,1]\to\mathbb{R}$ given by $g(t) = f(\gamma(t))$ to get the desired result:

$$f(b) - f(a) = g(1) - g(0) = g'(c) = \nabla f(\gamma(c))\cdot \gamma'(c) = 0$$

Also it's not stated that $f$ is differentiable, but since all partial derivatives $\nabla_i f = 0$ is continuous and $D$ is open this is guaranteed for all points in $D$ (see e.g. this question).

Solution 2:

Suppose $\nabla f=\left(\displaystyle\frac{\partial f}{\partial x},\displaystyle\frac{\partial f}{\partial y} \right)=(0,0)$.

$$ \frac{\partial f}{\partial x} = 0 \quad \Rightarrow \quad f(x,y)=A+g(y), $$ where $A\in \mathbb{R}$ and $g:\mathbb{R}\rightarrow \mathbb{R}$. Note that this only holds if $D$ is connected; e.g., if you consider $$f:]0,1[\cup ]1,2[\times \mathbb{R}\rightarrow \mathbb{R} \\(x,y)\mapsto\lfloor x \rfloor$$ it satisfies $\frac{\partial f}{\partial x} = 0$ but is not constant.

Then, we have $\frac{\partial f}{\partial y} =0= g'(y)$, therefore $ g(y)=B \in \mathbb{R} $. In summary $$ f(x,y)=A+B\in \mathbb{R} $$