Harmonic functions $u$ and $v$ such that $\nabla u = a \nabla v$ for some function $a$.

Since $a^2=(u_x^2+u_y^2)/(v_x^2+v_y^2)$, the function $a^2$ is differentiable on $\Omega$. Let $B$ be any open sub-rectangular region of $\Omega$ on which $a$ does not change its sign. Hence $a$ is differentiable on $B$. We claim that $a$ is a constant on $B$.

Differentiate $u_x=av_x$ with respect to $x$ to get $u_{xx}=a_xv_x+av_{xx}$. Similarly, $u_{yy}=a_yv_y+av_{yy}$. By comparing the two equations, we have $$a_xv_x=-a_yv_y.~~~~(1)$$ By differentiating $u_x=av_x$ with respect to $y$, we also $u_{xy}=a_yv_x+av_{xy}$, and similarly $u_{yx}=a_xv_y+av_{yx}$. Therefore, $$a_xv_y=a_yv_x.~~~~~~(2)$$ Square equations (1) and (2) and add. Since $v_x^2+v_y^2 \neq 0$, one has $a_x^2=a_y^2$. If $a$ is not constant on $B$, then $|a_x|=|a_y|\neq 0$ at some point. If at that point $a_x=a_y$, then $v_x=-v_y$ and $v_y=v_x$ hence $v_x=v_y=0$ a contradiction. Similarly, $a_x=-a_y$ leads to a contradiction. Hence $a$ is constant on $B$. By integrating $u_x=av_x$ on $B$, it follows that $u(x,y)=bv(x,y)+c$ on $B$ for some constants $b,c$. On any connected subset of $\Omega$ the constants $b$ and $c$ do not depend on $B$.

It seems that this works in any dimension.

Indeed, assume that $u,v$ are harmonic in a domain $\Omega\subseteq\mathbb R^n$, with $$\nabla u=a\nabla v$$ for a function $a$, and $\nabla v\neq 0$ everywhere in $\Omega$. Multiplying with $\nabla v$, we obtain that $a=\nabla u\nabla v|\nabla v|^{-2}$, which implies that $$|\nabla v|^2\nabla u=(\nabla u\nabla v)\nabla v.$$ Considering the divergence and using that $u,v$ are harmonic, we get $$\nabla\left(|\nabla v|^2\right)\nabla u=\nabla\left(\nabla u\nabla v\right)\nabla v.$$ On the other hand, the gradient of $a$ is equal to $$\nabla a=\frac{\nabla\left(\nabla u\nabla v\right)|\nabla v|^2-(\nabla u\nabla v)\nabla\left(|\nabla v|^2\right)}{|\nabla v|^4}=0,$$ so $a$ is constant.