I don't know the answers to your questions, i.e., I don't know whether there has been any work already done on this problem, nor do I know whether there is any kind of 'heat flow' argument for constructing solutions.

However, I suspect that the latter, if possible, is not going to be completely straightforward. After all, in the example I mention in my comment above, that of the $3$-sphere with $g_1$ being the standard metric with constant sectional curvature $1$ and $g_2=\lambda g_1$ for some constant $\lambda >1$, the surfaces of interest are the flat surfaces, of which there are many in the $3$-sphere, even compact ones near the Clifford torus. Moreover, because this is a hyperbolic problem in this case, one doesn't have good regularity for the 'stationary' solutions, so finding them by 'heat equation' methods appears to be a doubtful proposition.

In the general case (I'll assume that $M$ is oriented, for simplicity), one can understand the nature of the equations by using the structure equations on the oriented orthonormal frame bundle: Let $\pi:F\to M$ be the orthonormal frame bundle with respect to $g_1$. One has the tautological forms $\omega_i$ such that $\pi^*g_1 = {\omega_1}^2+{\omega_2}^2+{\omega_3}^2$ and the corresponding connection forms $\omega_{ij}=-\omega_{ji}$ satisfying the first structure equations
$$
d\omega_i = -\omega_{ij}\wedge \omega_j
$$
and the second structure equations
$$
d\omega_{ij} = -\omega_{ik}\wedge \omega_{kj} + \Omega_{ij}
=-\omega_{ik}\wedge \omega_{kj} + \tfrac12R_{ijkl}\ \omega_k\wedge\omega_l\ .
$$

Now, $F$ has a map to the unit sphere bundle $\nu:F\to\Sigma(M)$ defined by sending an orthonormal frame $(p;e_i)$ to the unit vector $(p,e_3)$. The $1$-form $\omega_3$ is the $\nu$-pullback of a well-defined $1$-form on $\Sigma(M)$, which, by abuse of notation, I write as $\omega_3$. Similarly, $\omega_1\wedge\omega_2$ is the $\nu$-pullback of a well-defined $2$-form on $\Sigma(M)$, which I denote by the same symbol. Meanwhile, although $\omega_{12}$ is not the $\nu$-pullback of any $1$-form on $\Sigma(M)$, its exterior derivative $d\omega_{12}= \omega_{31}\wedge\omega_{32} + \Omega_{12}$ *is* the $\nu$-pullback of a well-defined $2$-form on $\Sigma(M)$.

Now, if $S\subset M$ is any oriented surface in $M^3$, it has a Gauss map $\gamma:S\to \Sigma(M)$ given by sending a point of $S$ to its oriented normal vector. This map satisfies $\gamma^*\omega_3 = 0$, $\gamma^*(\omega_1\wedge\omega_2) = dA$, and $\gamma^*(d\omega_{12})=K\ dA$.

Now, let $g_2 = \lambda g_1$ and set $\lambda = e^{2\mu}$ for some $\mu$. Then the orthonormal frame bundle for $g_2$ can be compared with the orthonormal frame bundle of $g_1$ in the obvious way, and, under that identification, one has
$$
\omega^\ast_i = e^\mu\ \omega_i,\qquad\text{and}\qquad
\omega^\ast_{ij} = \omega_{ij} +\mu_j\omega_i - \mu_i\omega_j\ ,
$$
where I have decorated the forms associated to $g_2$ with an asterisk. In particular, one has, computing modulo $\omega_3$ (which is the same as computing modulo $\omega^\ast_3$),
$$
d\omega^\ast_{12} \equiv \omega_{31}\wedge\omega_{32}
-\mu_3(\omega_{31}\wedge\omega_2 + \omega_1\wedge\omega_{32})
+ (e^{2\mu}R^\ast_{1212}+{\mu_3}^2)\ \omega_1\wedge\omega_2\ .
$$
The condition that the surface $S$ have the same Gauss curvature with respect to $g_1$ and $g_2$ is then expressed as the condition that its Gauss map $\gamma$ pull back the $2$-form $\Upsilon = d\omega^\ast_{12}-e^{2\mu}d\omega_{12}$ to vanish, since $\gamma^\ast\Upsilon=K^\ast\ dA^\ast - e^{2\mu}K\ dA = e^{2\mu}(K^\ast-K)\ dA$. Now, one computes, using the above formulae, that
$$
\Upsilon = (1{-}e^{2\mu})\omega_{31}\wedge\omega_{32}
-\mu_3(\omega_{31}\wedge\omega_2{+}\omega_1\wedge\omega_{32})
+ \bigl(e^{2\mu}(R^\ast_{1212}{-}R_{1212}){+}{\mu_3}^2\bigr)\omega_1\wedge\omega_2\ .
$$

The exterior differential system generated by $\omega_3$, $d\omega_3$ and $\Upsilon$ is a Monge-Ampère system on the $5$-manifold $\Sigma(M)$. Its type depends on the sign of the quantity
$$
\Delta = (1-e^{2\mu})(R^\ast_{1212}-R_{1212}) - {\mu_3}^2.
$$
Where $\Delta>0$ it is hyperbolic, where $\Delta<0$ it is elliptic, and where $\Delta$ vanishes, it is degenerate.

**Note about P. Roitman's example:** It is interesting to note that this example is *not* elliptic everywhere. In fact, using the usual coordinates $x^1,x^2,x^3$ on $M=\mathbb{R}^3$ and letting $u = (u^1,u^2,u^3):\Sigma(M)\to S^2$ be the projection onto the unit $2$-sphere,
one has the formulae for the metrics he mentions in the form
$$
g_1 = (dx^1)^2+(dx^2)^2+(dx^3)^2\qquad\text{and}\qquad
g_2 = \frac{(dx^1)^2+(dx^2)^2+(dx^3)^2}{(x^3)^2}.
$$
Then the above formula for $\Delta$ works out to be
$$
\Delta = \frac{1-(x^3)^2-(u^3)^2}{(x^3)^2}.
$$
In particular, the equation is elliptic in the entire region $x^3 >1$, but in the slab $0 < x^3 < 1$, only the solution surfaces whose normals are 'close enough' to vertical are defined by an elliptic equation. In other words, there are hyperbolic solutions in this 'boundary region'. It might be interesting to see what these non-analytic solutions look like.

*Added Comment 1:* This example has degenerate solutions (i.e., solutions whose $1$-graphs lie in the locus $\Delta=0$) depending on an arbitrary function of one variable: In fact, one can parametrize them in the form
$$
x^1=a(s)+\cos(s)\bigl(t{-}\tanh t \bigr), \quad
x^2=b(s)+\sin(s)\bigl(t{-}\tanh t \bigr), \quad
x^3 = \text{sech}\ t
$$
where $a$ and $b$ are functions of $s$ that satisfy $a'(s)\cos(s)+b'(s)\sin(s)=0$.
Choosing $a$ and $b$ to be non-analytic produces non-analytic surfaces.

*Added Comment 2:* It turns out that there is a Lagrangian $\Lambda_E$ for surfaces of elliptic type (i.e., ones for which $\Delta<0$) in the upper half space such that the elliptic integral surfaces of the above system are extrema of $\Lambda_E$, and there is a Lagrangian $\Lambda_H$ for surfaces of hyperbolic type (i.e., ones for which $\Delta>0$) in the upper half space such that the hyperbolic integral surfaces of the above system are extrema of that $\Lambda_H$. Each of these Lagrangians blows up along the locus $\Delta=0$ in $\Sigma(M)$, and it appears that the degenerate solutions listed in the Comment 1 above are not the extrema of any Lagrangian.