By property 1, the function $\psi_n$ increases going inside the domain. You should set
$$\nu(x_0) = - \frac{\nabla \psi_n(x_0)}{||\nabla \psi_n(x_0)||}$$
if you want the outward normal.
How do we know that $\nabla \psi_n (x_0) \neq 0$?
By assumption, $\psi$ is a diffeomorphism (its inverse is also smooth). Therefore, the Jacobian matrix of $\psi$ is invertible. Consequently, every row of this matrix is a nonzero vector. Including the last row, which is $\nabla \psi_n$.
how do we know that $\nu$ is well defined?
Suppose $\phi$ is another map like $\psi$. It's convenient to normalize it (by adding a constant) so that $\phi(x_0)=\psi(x_0)$. Then $\phi = F\circ \psi$ where $F = \phi\circ \psi^{-1}$ is a diffeomorphism of a neighborhood of $\psi(x_0)$. Note that $F$ leaves the hyperplane $x_n=0$ invariant. Therefore, the partial derivatives of its $n$th component with respect to $x_1,\dots,x_{n-1}$ are zero. In matrix terms, this means the $n$th row of the Jacobian matrix $DF$ is of the form $(0,0,\dots, 0, * )$ where $*$ stands for a nonzero number. This number is positive, because $F$ sends upper halfspace into itself.
The chain rule, applied to $\phi = F\circ \psi$, shows that $\nabla \phi_n(x_0)$ is a positive multiple of $\nabla \psi_n(x_0)$. Therefore, the normalization of these vectors produces the same result.
The conditions under which the divergence theorem holds can be much weaker than what you propose. We do not even need that the partial derivatives of the vector field be continuous on $\Omega$ -- bounded and Lebesgue integrable is sufficient.
There are even weaker forms. This paper by Bochner is a good starting point as a reference.
To see that continuity of the partial derivatives on the boundary is not necessary, consider the example (for the analogous Green's theorem) where $\Omega = [0,1]^2$ and
$$P(x,y) = 0 \\ Q(x,y) = \begin{cases}yx^2\sin(1/x), \, \, x \neq 0 \\ 0, \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, x = 0 \end{cases}$$
In this case, we have
$$\iint_{\Omega} \left(\frac{\partial Q}{\partial x}-\frac{\partial P}{\partial y} \right)dA = \int_{\partial \Omega} P \,dx+Q\,dy $$
even though $\frac{\partial Q}{\partial x}$ is not continuous at points in $\{(x,y): x = 0, 0 < y \leqslant 1 \}$.
Best Answer
As suggested by fourierwho in their comment, perhaps the most the natural domains for which the divergence (also called Gauss-Green) theorem holds are the sets of finite perimeter, i.e. Caccioppoli sets, so let's precisely see why.
Definition 1 ([1], §3.3 p. 143). Let $\Omega$ a Lebesgue measurable set in $\mathbb{R}^n$. For any open subset $G\subseteq\mathbb{R}^n$ the perimeter of $\Omega$ in $G$, denoted as $P(\Omega,G)$, is the variation of $\chi_\Omega$ in $\Omega$ i.e. $$ \begin{split} P(\Omega,G)&=\sup\left\{\int_\Omega \nabla\cdot\varphi\,\mathrm{d}x\,:\,\varphi\in [C_c^1(G)]^n, \|\varphi\|_\infty\leq1\right\}\\ & =| \nabla \chi_{\Omega\cap G}|=TV(\Omega,G) \end{split}\tag{1}\label{1} $$ where $[C_c^1(G)]^n$ is the set of compact support continuously differentiable vector functions in $G$ and $TV$ is the total variation of the set function $\nabla \chi_{\Omega\cap G}$.
The set $\Omega$ is a set of finite perimeter (a Caccioppoli set) in $G\subseteq\mathbb{R}^n$ if $P(\Omega,G)<\infty$.
Why definition \eqref{1} implies a natural extension of the classical divergence (Gauss-Green) theorem? For simplicity lets consider sets of finite perimeter: $P(\Omega)<\infty$ implies that the distributional derivative of the characteristic function of $\Omega$ is a vector Radon measure whose total variation is the perimeter defined by \eqref{1}, i.e. $$ \nabla\chi_\Omega(\varphi)=\int_\Omega\nabla\cdot\varphi\,\mathrm{d}x=\int_\Omega \varphi\,\mathrm{d}\nabla\chi_\Omega\quad \varphi\in [C_c^1(\mathbb{R}^n)]^n\tag{2}\label{2} $$ Now the support in the sense of distributions of $\nabla\chi_\Omega$ is $\subseteq\partial\Omega$ ([2], §1.8 pp. 6-7): to see this note that if $x\notin\partial\Omega$, it should belong to an open set $A\Subset\mathbb{R}^n$ such that either $A\Subset\Omega$ or $A\Subset\mathbb{R}^n\setminus\Omega$:
Also, as a general corollary of (one of the versions of) Radon-Nikodym theorem ([1], §1.1 p. 14) we can apply a polar decomposition to $\nabla\chi_\Omega$ and obtain $$ \nabla\chi_\Omega=\nu_\Omega|\nabla\chi_\Omega|_{TV}\equiv\nu_\Omega|\nabla\chi_\Omega|\tag{3}\label{3} $$ where $\nu_\Omega$ is a $L^1$ function taking values on the unit sphere $\mathbf{S}^{n-1}\Subset\mathbb{R}^n$, and rewriting \eqref{2} by using \eqref{3} we obtain the sought for general divergence (Gauss-Green) theorem $$ \int_\Omega\!\nabla\cdot \varphi\, \mathrm{d}x =\int_{\partial\Omega} \!\varphi\,\cdot\nu_\Omega\, \mathrm{d}|\nabla\chi_\Omega|\quad\forall\varphi\in [C_c^1(\mathbb{R}^n)]^n\tag{4}\label{4} $$ Note that this result is an almost direct consequence of definition 1 above, with minimal differentiability requirement imposed on the data $\varphi$: it seems to follow directly from the given definition of perimeter \eqref{2} through the application of general (apparently unrelated) theorems on the structure of measures and distributions, and in this sense it is the most "natural form" of the divergence/Gauss-Green theorem.
Further notes
[1] Ambrosio, Luigi; Fusco, Nicola; Pallara, Diego (2000), Functions of bounded variation and free discontinuity problems. Oxford Mathematical Monographs, New York and Oxford: The Clarendon Press/Oxford University Press, New York, pp. xviii+434, ISBN 0-19-850245-1, MR1857292, Zbl 0957.49001.
[2] Giusti, Enrico (1984), Minimal surfaces and functions of bounded variations, Monographs in Mathematics, 80, Basel–Boston–Stuttgart: Birkhäuser Verlag, pp. XII+240, ISBN 978-0-8176-3153-6, MR 0775682, Zbl 0545.49018