I think that the answer is 'yes' if $U$ is simply-connected, because there is a way to construct a candidate 'approximate surface' from 'approximate solutions' of Gauss and Codazzi, but a more useful answer would be one that gave you estimates of how close the metric and second fundamental form of the approximate surface are to the given data in terms of the failure of the Gauss and Codazzi equations of the original data.
Here is what I mean by constructing an approximate solution: There's no loss of generality in assuming that $U$ is the set in the $uv$-plane given by $|u|,|v|\le 1$ and that $g$ and $l$ are continuous or even smooth on the boundary. We can choose $1$-forms $\omega_1$ and $\omega_2$ on $U$ such that $g = \omega_1^2+\omega_2^2$ and we can write $l = l_{ij}\omega_i\omega_j$ where $l_{ij}=l_{ji}$ and set $\omega_{3i} = l_{ij}\omega_j$. Also, let $\omega_{12}=-\omega_{21}$ satisfy $\mathrm{d}\omega_i = -\omega_{ij}\wedge\omega_j$.
Now consider the matrix valued $2$-form
$$
\eta = \begin{pmatrix}
0 & 0 & 0 & 0\\
\omega_1 & 0 &\omega_{12} & -\omega_{31}\\
\omega_2 & -\omega_{12} & 0 & -\omega_{32}\\
0 & \omega_{31} & \omega_{32} & 0
\end{pmatrix} = A\,\mathrm{d}u + B\,\mathrm{d} v.
$$
The Gauss and Codazzi equations would hold if and only if
$\mathrm{d}\eta + \eta\wedge\eta = 0$, i.e.
$$
B_u-A_v + [A,B]=0.
$$
In general, we'll have
$$
\mathrm{d}\eta + \eta\wedge\eta
=\begin{pmatrix}
0 & 0 & 0 & 0\\
0 & 0 &H & -C_1\\
0 & -H & 0 & -C_2\\
0 & C_1 & C_2 & 0
\end{pmatrix}\,\mathrm{d}u\wedge\mathrm{d}v
= F\,\mathrm{d}u\wedge\mathrm{d}v \,,
$$
where $C_1$, $C_2$, and $H$ are 'small' in some appropriate sense (which it will be important to make precise if you want this to be useful).
Now, consider the solution of the matrix ODE with initial condition
$$
L'(u) = L(u)A(u,0), \qquad L(0)=I_4\,.
$$
and use this to solve the matrix ODE with initial condition
$$
M_v(u,v) = M(u,v)B(u,v),\qquad M(u,0) = L(u).
$$
Of course, $M$ won't satisfy $M^{-1}\mathrm{d}M = \eta$ unless $C_1=C_2=H=0$.
But, by construction, $M^{-1}\mathrm{d}M-\eta = E\,\mathrm{d}u$, where it should be possible to estimate the size of $E$ in terms of the given data $g$, $l$ and $C_1$, $C_2$, $H$, and this estimate should be boundable in such a way that when $C_1$, $C_2$, and $H$ are small in the appropriate sense, then $E$ will be as well. The point is that $E(u,v)$
satisfies the matrix ODE with initial condition
$$
E_v = [E,B] + F,\qquad E(u,0) = 0,
$$
so $|E|$ can be bounded in terms of $|B|$ and $|F|$. (Clearly, if $F\equiv0$, then $E\equiv0$.)
The first column of $M$ will look like
$$
\begin{pmatrix}1\\x(u,v)\\y(u,v)\\z(u,v)\end{pmatrix},
$$
and $(x,y,z):U\to\mathbb{R}^3$ will give an 'approximate surface'. With the right estimates on $C_1$, $C_2$, and $H$ (which may involve coefficients that depend on $g$ and $l$), you should be able to say how close its first and second fundamental form are to the originally given data.
Of course, obtaining good estimates would be an interesting project. What one would really want would be estimates that were expressed in terms of the $g$-geometry of the open set $U$ and the $g$-geometry of the second fundamental form, without reference to the local coordinates $(u,v)$. This should be possible, and I think that the estimates would be interesting. I suspect something like this has been done in the literature, maybe in the computer science literature where it might be used in surface reconstruction from local data (which is only approximately known). I wouldn't be surprised to find that there is something about this in the discrete differential geometry literature, but I'm not familiar with much of that enormous corpus.