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.
