Let $U$ be an open neighbuorhood $U \subset \mathbb{R}^2$, and let $g$ be a smooth Riemannian metric on $U$.
Let $x$ be a smooth function on $U$, with nonvanishing derivative $dx \neq 0$.
Let $p \in U$. Does there always exist a smooth function $y$ on a neighbourhood of $p$, such that the coordinates $x,y$ are orthogonal w.r.t the metric $g$?
i.e. $g=f(x,y)dx^2+h(x,y)dy^2$, for some functions $f,h$.
Equivalently: $g(\partial_x,\partial_y)=0$.
Note that I am not asking for Isothermal coordinates, so $f,h$ can be different.
Yes, there always exists such a function in a neighborhood of each point.
Given $x$, let $\operatorname{grad}_g x$ denote the gradient of $x$ with respect to $g$. Another function $y$ makes an orthogonal coordinate system $(x,y)$ on some neighborhood of $p\in U$ if and only if $dy_p\ne 0$ and $\operatorname{grad}_g x$ is in the kernel of $dy$ at each point (exercise).
The canonical form theorem for nonvanishing vector fields (also called the straightening theorem) shows that there are coordinates $(u,v)$ on a neighborhood of $p$ such that $\operatorname{grad}_g x \equiv \partial/\partial u$ on that neighborhood. In particular, this means that $\operatorname{grad}_g x$ is in the kernel of the $1$-form $dv$ everywhere on that neighborhood, so $y=v$ is the desired coordinate function.