# Constitutive relationships

## Water retention curve, van Genuchten (1980):

\begin{aligned} \theta(\psi) = \theta_{r} + (\theta_{s} - \theta_{r}) (1+(a_{VG} \psi)^{n_{VG}})^{-m_{VG}} \\\end{aligned}

## Hydraulic conductivity function, van Genuchten et al. (1991):

\begin{aligned} k(\psi) = k_{sat} \frac {(1-((a_{VG} \psi)^{n_{VG}m_{VG}}) (1+(a_{VG} \psi)^{n_{VG}})^{-m_{VG}}))^2} { (1+(a_{VG} \psi)^{n_{VG}})^{m_{VG}l_{VG}}}\end{aligned}

# Partial differential equation

Mass conservation law.

\begin{aligned} \frac{\partial \theta}{\partial t} = - \nabla \cdot \vec{q} \\\end{aligned}

Darcy law.

\begin{aligned} \vec{q} = k \left( h, P \right) \nabla H \\\end{aligned}

Where h is the pore pressure, P is the position in 3D space and H is the total head.

\begin{aligned} P = (x, y, z)\\ h= - \psi \\ \nabla H = \nabla h + \nabla z \\\end{aligned}

The partial differential equation can be written as:

\begin{aligned} \frac{\partial \theta}{\partial t} = -\nabla \cdot k \left( h, P \right) \nabla H + A, \\\end{aligned}

where $$A$$ is a source/sink term.

# Weak formulation

Integrate on both sides and multiply by the test function $$v$$.

\begin{aligned} \int_{\Omega}\frac{\partial \theta}{\partial t} v\,d\Omega = -\int_{\Omega}\nabla \cdot k \left( h, P \right) \nabla H v\,d\Omega + \int_{\Omega}A v\,d\Omega \\\end{aligned}

Integrate by parts.

\begin{aligned} \int_{\Omega}\nabla \cdot k \left( h, P \right) \nabla H v\,d\Omega = -\int_{\Omega} k \left( h, P \right) \nabla H \nabla v\,d\Omega + \\ \int_{\Omega}\nabla \cdot \left( k \left( h, P \right) \nabla H v\right) \,d\Omega \\\end{aligned}

Apply Gauss divergence theorem.

\begin{aligned} \int_{\Omega}\nabla \cdot k \left( h, P \right) \nabla H v\,d\Omega = -\int_{\Omega} k \left( h, P \right) \nabla H \nabla v\,d\Omega + \int_{\Gamma_{N} \cup \Gamma_{D}} k \left( h, P \right) \nabla H v n\,d\Gamma \\ \int_{\Omega}\nabla \cdot k \left( h, P \right) \nabla H v\,d\Omega = -\int_{\Omega} k \left( h, P \right) \nabla H \nabla v\,d\Omega + \int_{\Gamma_{N}} k \left( h, P \right) \nabla H v n\,d\Gamma + \int_{\Gamma_{D}} k \left( h, P \right) \nabla H v n\,d\Gamma\end{aligned}

Because $$v = 0$$ on $$\Gamma_{D}$$.

\begin{aligned} \int_{\Omega}\nabla \cdot k \left( h, P \right) \nabla H v\,d\Omega = -\int_{\Omega} k \left( h, P \right) \nabla H \nabla v\,d\Omega + \int_{\Gamma_{N}} k \left( h, P \right) \nabla H v n\,d\Gamma \\\end{aligned}

Put in the initial integral.

\begin{aligned} \int_{\Omega} \frac{\partial \theta}{\partial t} v\,d\Omega = \int_{\Omega} k \left( h, P \right) \nabla H \nabla v\,d\Omega - \int_{\Gamma_{N}} \left( k \left( h, P \right) \nabla H n \right) v\,d\Gamma + \int_{\Omega}A v\,d\Omega \\\end{aligned}

$$\theta \left( h \right)$$ is a constitutive relationship. The weak formulation can thus be written on a $$h$$ basis.

\begin{aligned} \int_{\Omega} \frac{\partial \theta}{\partial t} \times \frac{\partial h}{\partial h} v\,d\Omega = \int_{\Omega} k \left( h, P \right) \nabla H \nabla v\,d\Omega - \int_{\Gamma_{N}} \left( k \left( h, P \right) \nabla H n \right) v\,d\Gamma + \int_{\Omega}A v\,d\Omega \\ \\ \frac{\partial \theta}{\partial h} \int_{\Omega} \frac{\partial h}{\partial t} v\,d\Omega = \int_{\Omega} k \left( h, P \right) \nabla H \nabla v\,d\Omega - \int_{\Gamma_{N}} \left( k \left( h, P \right) \nabla H n \right) v\,d\Gamma + \int_{\Omega}A v\,d\Omega \\\end{aligned}

# Notes

Integration by parts.

\begin{aligned} \int_{\Omega} \left(\nabla \cdot u \right) v\,d\Omega = -\int_{\Omega} u \nabla v\,d\Omega + \int_{\Omega} \nabla \cdot \left( u v \right) n\,d\Omega \\\end{aligned}

Gauss divergence theorem.

\begin{aligned} \int_{\Omega} \nabla \cdot F\, d\Omega = \int_{\partial \Omega = \Gamma_{N} \cup \Gamma_{D}} F \cdot n\, d\Gamma\end{aligned}

The problem in steady-state with no Neumann boundary conditions can be defined in Sfepy.