User:Jdavi231/sandbox

Solution Methods
Due to the lack of analytical closed form solutions to Richard's equation, numerical methods must be used to solve the equation. A commonly used numerical method is the finite difference method. In this case, we will show the finite difference formulation of the mixed form of Richards Equation, using backward Euler method for the time discretization. In general, solutions using the mixed form of Richard's equation result in a better mass conservation than other forms. Due to the non linear-behavior of Richards equation, iterative solutions methods such as Picard iteration are utilized.

Beginning with the 1D mixed form of Richards equation we have

$$\frac{\partial \theta}{\partial t}=\frac{\partial}{\partial z} \left[K \left ( \frac{\partial h}{\partial z} \right) \right] +\frac{\partial K }{\partial z} $$

Applying a backward Euler approximation for the time discretization

$$\frac{\partial \theta}{\partial t} \approx \frac{\theta_i^{n+1,m+1} - \theta_i^{n}}{\Delta t}$$

where $$\Delta t $$ denotes the size of the time step, $$n $$ denotes the time step ($$t = n\Delta t $$), m denotes the iteration level of the solution at that time step, and index $$i $$ denotes the index of the solution node ($$z = i\Delta z $$). Applying a finite difference approximation to the remaining terms on the right hand side, we have

$$\frac{\partial K}{\partial z} \approx \frac{K_{i+1/2}^{n+1,m} - K_{i-1/2}^{n+1,m} }{\Delta z}$$

$$\frac{\partial}{\partial z} \left[K \left ( \frac{\partial h}{\partial z} \right) \right] \approx \frac{K_{i+1/2}^{n+1,m}}{(\Delta z)^2}(h_{i-1}^{n+1,m+1}-h_i^{n+1,m+1}) + \frac{K_{i-1/2}^{n+1,m}}{(\Delta z)^2}(h_{i+1}^{n+1,m+1}-h_i^{n+1,m+1}) $$

where $$\Delta z $$ denotes the spacing between solution nodes and $$K_{i\pm1/2} $$is the hydraulic conductivity at the interface between nodes. This is usually taken to be the geometric mean of the Hydraulic conductivity calculated at node $$i $$ and node $$i\pm1 $$. We are now left with the following equation

$$ \frac{\theta_i^{n+1,m+1} - \theta_i^{n}}{\Delta t} = \frac{K_{i+1/2}^{n+1,m}}{(\Delta z)^2}(h_{i-1}^{n+1,m+1}-h_i^{n+1,m+1}) + \frac{K_{i-1/2}^{n+1,m}}{(\Delta z)^2}(h_{i+1}^{n+1,m+1}-h_i^{n+1,m+1}) + \frac{K_{i+1/2}^{n+1,m} - K_{i-1/2}^{n+1,m} }{\Delta z} $$

In order to cast this equation into a readily solvable form, we approximate $$\theta_i^{n+1,m+1}$$ to first order about $$h_i^{n+1,m}$$.

$$\theta_i^{n+1,m+1} \approx \theta_i^{n+1,m} + \frac{\partial\theta}{\partial h}|_{n+1,m} (h_i^{n+1,m+1}-h_i^{n+1,m}) + O(\delta^2)$$

In most literature the derivative $$\frac{\partial\theta}{\partial h}|_{n+1,m}$$is written as $$C_i^{n+1,m}$$. Finally, applying the substitution $$\delta^{n+1,m+1}_i = h^{n+1,m+1}_i - h^{n+1,m}_i$$we are left with

$$C_i^{n+1,m}\frac{\delta_i^{n+1,m}}{\Delta t} + \frac{K_{i-1/2}^{n+1,m}}{(\Delta z)^2}(\delta_{i-1}^{n+1,m+1} - \delta_i^{n+1,m+1}) + \frac{K_{i+1/2}^{n+1,m}}{(\Delta z)^2}(\delta_{i+1}^{n+1,m+1} - \delta_i^{n+1,m+1}) = $$

$$\frac{K_{i-1/2}^{n+1,m}}{(\Delta z)^2}(h_{i-1}^{n+1,m} - h_i^{n+1,m})+ \frac{K_{i+1/2}^{n+1,m}}{(\Delta z)^2}(h_{i+1}^{n+1,m} - h_i^{n+1,m}) - \frac{K_{i+1/2}^{n+1,m} - K_{i-1/2}^{n+1,m}}{(\Delta z)^2} -\frac{\theta_i^{n+1,m}- \theta_i^n}{\Delta t}$$

This system of equations can be written as a matrix equation of the form $$Ax = b $$ and solved for x which is a vector of all $$\delta_i^{n+1,m+1}$$. The solution of this system of nonlinear equations is computationally simple, due to matrix A being tridiagonal. In the case of a tridiagonal matrix the Thomas algorithm can be implemented to solve for x. In practice this system is solved iteratively for all $$\delta_i^{n+1,m+1}$$until the magnitude of all the $$\delta$$ are less than some chosen tolerance.

Richards Equation Coupled to Heat Transfer
In some instances it may not be justified to treat flow in porous media as being purely due to capillary action. Water transfer in unsaturated soil is also caused by temperature gradients in the material. In this case the only gradients to consider are temperature gradients in the vertical direction since soil with horizontal temperature gradients are rarely seen. Starting with the 1D form of the conservation of mass equation from earlier we have

$$\frac{\partial \theta}{\partial t}= -\frac{\partial}{\partial z} q$$

where q is the associated flux in the z direction. Including flow due to temperature gradients as well the Darcy flux

$$q= q_{Darcy} + q_{Temp} = - K \frac{\partial H}{\partial z} - D_T\frac{\partial T}{\partial z}$$

substituting this new q into the mass conservation equation we get:

$$\frac{\partial \theta}{\partial t} = \frac{\partial}{\partial z} [K \frac{\partial h}{\partial z}] + \frac{\partial}{\partial z}[D_T\frac{\partial T}{\partial z}] + \frac{\partial K}{\partial z}$$

where $$D_T$$ is the moisture diffusivity under the influence of a temperature gradient. We now have a coupled equation for unsaturated flow that includes effects from capillary action and temperature gradients.

Helmholtz Free Energy
Under the assumption that the soil is in an open thermodynamic system, the Helmholtz free energy is an appropriate thermodynamic potential for describing the system. The differential change in the Helmholtz Free energy $$F$$ is therefore given as

$$dF = -SdT - \rho dV + \mu dM$$

where S is the entropy, T is the temperature of the system, $$\rho$$ is the pressure, $$V $$ is the soil volume, and $$\mu$$ is the chemical potential. The water mass $$M$$ = $$\rho \theta V$$ where $$\rho$$ is the density of water. Using Richards equation to describe movement of water in the soil, we can assume that any soil water fluxes are isotherm processes and that no expansion, or settling of the soil occurs. With this assumption we can take $$SdT$$ and $$\rho dV$$ to be zero. We can then calculate the Helmholtz free energy as

$$F = V\rho \int\mu d\theta + c$$

substituting the $$\mu$$ with pressure head $$H$$= $$h + z$$ we get:

$$F = \rho V \int h(\theta)d \theta + \int z d\theta +c $$