I am interested in the following three surfaces in $\mathbb{R}^3$:
$$S_1=\{(x,y,x+y+x^3+\sqrt{y}): x \in [0,x_0], y\in [0,y_0]\},$$
$$S_2=\{(x,y,x+y+x^2 +y^2+1): x \in [0,x_0], y\in [0,y_0]\},$$
$$S_3=\{(x,y,x+y+x^4 +y^3): x \in [0,x_0], y\in [0,y_0]\}$$ for some $x_0, y_0>0$.
I wonder if there is a natural extension of the notion of Lebesgue measure (in $\mathbb{R}^n$) on these surfaces and a natural notion of open sets on these surfaces such that something like the Lebesgue density theorem would hold using these two extended notions?
What you want is the area formula from multivariable calculus. In your case each of the three surfaces $S_1,S_2,S_3$ is the graph of a function $z=f(x,y)$ defined for all $(x,y)$ in some (measurable) subset $D$ of the plane. Any (measurable) subset of the surface can be written in the form $f(A)$ for some (measurable) subset $A \subset D$. Furthermore, $f(A)$ is an open subset of the surface if and only if $A$ is an open subset of $D$. The area formula gives us $$Area(f(A)) = \int\!\!\int_A \sqrt{\left(\frac{\partial f}{\partial x}\right)^2 + \left(\frac{\partial f}{\partial y}\right)^2 + 1} \,\,\,\, dx \, dy $$ and that's what you can use as Lebesgue measure on the surface. The Lebesgue density theorem will hold for this measure.