The problem is that, in general, the faces of an irregular hexahedron are not flat.
There is not one normal, the surface area is curved in space and difficult to
calculate exactly (whatever the additional value of "exact" is in a numerical context).
So I think you should anyway be satisfied with some sensible approximations. Here comes a proposal.
Let the vertices of the (in general non-planar) face of an irregular hexahedron be given
by $\,\vec{r}_A , \vec{r}_B , \vec{r}_C , \vec{r}_D$ ,
as shown in the picture. Then form the vertices of the gray (
Varignon) quadrilateral as :
$$\vec{r}_1 = \frac{1}{2}\left(\vec{r}_A+\vec{r}_C\right) \quad ; \quad
\vec{r}_2 = \frac{1}{2}\left(\vec{r}_B+\vec{r}_D\right) \\
\vec{r}_3 = \frac{1}{2}\left(\vec{r}_A+\vec{r}_B\right) \quad ; \quad
\vec{r}_4 = \frac{1}{2}\left(\vec{r}_C+\vec{r}_D\right)
$$
It can be easily proved that the quadrilateral $\overline{1324}$ is a paralellogram and hence
it is planar. The intersection point of the (yellow) diagonals of the Varignon paralellogram
is a good place to define the normal, as a vector perpendicular to $\overline{1324}$ . Last but not
least the total area of the face is defined by the sum of five planar areas: the area of the
paralellogram $\overline{1324}$ and the four areas of the triangles $\,\overline{A31} ,
\overline{B23} , \overline{D42} , \overline{C14}$ . Hope this helps.
Lets examine the quadrilateral as a linearly interpolated surface. If the vertices are
$$\bbox{
\vec{p}_1 = \left [ \begin{matrix} x_1 \\ y_1 \\ z_1 \end{matrix} \right ]
} , \quad \bbox {
\vec{p}_2 = \left [ \begin{matrix} x_2 \\ y_2 \\ z_2 \end{matrix} \right ]
} , \quad \bbox {
\vec{p}_3 = \left [ \begin{matrix} x_3 \\ y_3 \\ z_3 \end{matrix} \right ]
} , \quad \bbox {
\vec{p}_4 = \left [ \begin{matrix} x_4 \\ y_4 \\ z_4 \end{matrix} \right ]
}$$
and we parametrise the surface using a vector-valued function $\vec{p}(u, v)$, with $0 \le u \le 1$, $0 \le v \le 1$, then
$$\bbox{ \vec{p}(u, v) = (1 - v)\biggr( (1-u) \vec{p}_1 + u \vec{p}_2 \biggr) + v \biggr( (1-u) \vec{p}_3 + u \vec{p}_4 \biggr) }$$
i.e.
$$\bbox{ \vec{p}(u, v) = \vec{p}_1 + u \left( \vec{p}_2 - \vec{p}_1 \right) + v \left( \vec{p}_3 - \vec{p}_1 \right) + u v \left( \vec{p}_4 - \vec{p}_3 + \vec{p}_1 - \vec{p}_2 \right) }$$
If the four vertices are nonplanar, the surface is curved. In all cases, the surface passes through
$$\bbox{ \vec{p}(0.5, 0.5) = \frac{ \vec{p}_1 + \vec{p}_2 + \vec{p}_3 + \vec{p}_4 }{4} }$$
The surface tangent vectors are
$$\bbox{ \begin{aligned}
\vec{p}_u (u , v) &= \frac{\partial \vec{p}(u, v)}{\partial u} = \vec{p}_2 - \vec{p}_1 + v \left( \vec{p}_4 - \vec{p}_3 - \vec{p}_2 + \vec{p}_1 \right)
\\
\vec{p}_v (u, v) &= \frac{\partial \vec{p}(u, v)}{\partial v} = \vec{p}_3 - \vec{p}_1 + u \left( \vec{p}_4 - \vec{p}_3 - \vec{p}_2 + \vec{p}_1 \right) \\
\end{aligned} }$$
and the exact area of this surface is
$$\bbox{ A = \int_0^1 d u \int_0^1 d v \; \vec{p}_u (u, v) \cdot \vec{p}_v (u, v) = \frac{ \left( \vec{p}_1 + \vec{p}_2 - \vec{p}_3 - \vec{p}_4 \right)\cdot\left( \vec{p}_1 - \vec{p}_2 + \vec{p}_3 - \vec{p}_4 \right) }{4} }$$
The surface normal vector is $\vec{n}(u, v)$,
$$\bbox{ \vec{n}(u , v) = \frac{\partial \vec{p}(u, v)}{\partial u} \times \frac{\partial \vec{p}(u, v)}{\partial v} }$$
and the mean (expected value) of the normal vector is
$$\bbox{ \langle\vec{n}\rangle = \int_0^1 d u \int_0^1 d v \; \vec{n}(u ,v) = \left [ \begin{matrix}
\frac{(y_4 - y_1)(z_3 - z_2) - (z_4 - z_1)(y_3 - y_2)}{2} \\
\frac{(z_4 - z_1)(x_3 - x_2) - (x_4 - x_1)(z_3 - z_2)}{2} \\
\frac{(x_4 - x_1)(y_3 - y_2) - (y_4 - y_1)(x_3 - x_2)}{2} \\
\end{matrix} \right] = \vec{n}\left(\frac{1}{2} , \frac{1}{2}\right) }$$
If we consider a flow through the surface, we could assume the flow is perpendicular to the surface normal. Then, the effective surface is the area of the surface when projected to a plane perpendicular to the unit normal vector $\hat{n}$ used,
for example
$$\hat{n} = \frac{\langle\vec{n}\rangle}{\left\lVert\langle\vec{n}\rangle\right\rVert}$$
$$A_\perp = \int_0^1 d u \int_0^1 d v \; \Biggr( \vec{p}_u (u, v) - \hat{n} \biggr(\hat{n}\cdot\vec{p}_u (u,v) \biggr) \Biggr) . \Biggr( \vec{p}_v (u, v) - \hat{n} \biggr(\hat{n}\cdot\vec{p}_v (u,v) \biggr) \Biggr)$$
It turns out that for this choice of $\hat{n}$, $A_\perp = A$.
Therefore, choosing
$$\bbox[#ffffef, 1em]{ \vec{n} = \left [ \begin{matrix}
(y_4 - y_1)(z_3 - z_2) - (z_4 - z_1)(y_3 - y_2) \\
(z_4 - z_1)(x_3 - x_2) - (x_4 - x_1)(z_3 - z_2) \\
(x_4 - x_1)(y_3 - y_2) - (y_4 - y_1)(x_3 - x_2) \end{matrix} \right ] = \; \bigr( \vec{p}_4 - \vec{p}_1 \bigr) \times \bigr( \vec{p}_3 - \vec{p}_2 \bigr) }$$
$$\bbox{\hat{n} = \frac{\vec{n}}{\left\lVert\vec{n}\right\rVert} = \frac{\vec{n}}{\sqrt{\vec{n} \cdot \vec{n}}}}$$
$$\bbox[#ffffef, 1em]{ A = \frac{ \left( \vec{p}_1 + \vec{p}_2 - \vec{p}_3 - \vec{p}_4 \right)\cdot\left( \vec{p}_1 - \vec{p}_2 + \vec{p}_3 - \vec{p}_4 \right) }{4} }$$
is an obvious choice, in my opinion.
Best Answer
Take $\eta = 0$ for instance. A tangent vector is $$ \left[ \begin{matrix} \frac{\partial x}{\partial \xi}(\xi, 0) \\ \frac{\partial y}{\partial \xi}(\xi, 0) \end{matrix} \right] $$ By rotation of -90 degrees, an outward normal vector is $$ \left[ \begin{matrix} - \frac{\partial y}{\partial \xi}(\xi, 0) \\ \frac{\partial x}{\partial \xi}(\xi, 0) \end{matrix} \right] $$ You can proceed similarly for the other edges.