Issues with Richards Equation Simulation

Hello FEniCS community,

I am working on a simulation involving the Richards equation for unsaturated flow in soils with fenics 2019.1.0 installed from conda-forge. The general form of Richards’ equation can be written as:

\frac{\partial \theta(u)}{\partial t} = \nabla \cdot [K(u) (\nabla u + \nabla z)]


  • \theta(u) is the volumetric water content as a function of pressure head u.

  • K(u) is the hydraulic conductivity as a function of pressure head u.

  • z is the elevation.

For variational formulation,
\int_{\Omega} \frac{\partial \theta(u)}{\partial t} v ~ d\Omega - \int_{\Omega} K(u) \nabla (u + z) \cdot \nabla v ~ d\Omega = \int_{\partial\Omega} v K(u) \nabla (u + z) \cdot \mathbf{n} ~ ds


\int_{\Omega} \frac{\theta(u) - \theta(u_{old})}{\Delta t} v ~ d\Omega - \int_{\Omega} K(u) \nabla (u + z) \cdot \nabla v ~ d\Omega - \int_{\partial\Omega} v q ~ ds = 0

Here, u_{old} is the pressure head from the previous time step and \Delta t is the time step size.

My expectation is that the pressure head, u, would converge over time to a profile where u decreases linearly with elevation (x[1]), then breaks to a minimal value of u such as K(u) = q. My current results don’t match this expectation. Below is my complete code for the simulation:

import numpy as np
from fenics import *
import matplotlib.pyplot as plt

# Soil parameters
a = 0.2
theta_r = 0.05
theta_s = 0.45
n = 2.0
m = 1-1/n
k_sat = 1E-3
l = 0.5

# Soil water content and hydraulic conductivity
theta = lambda u: conditional(lt(u, 0), theta_r + (theta_s - theta_r) / (1 + (-a * u) ** n) ** m, theta_s)
K = lambda u: conditional(lt(u, 0), k_sat * (1 - (-a * u)**(n-1) * (1 + (-a * u)**n)**(-m)) * (1 + (-a * u)**n)**(-m*l), k_sat)

# Mesh and function space
ymin, ymax = 0, 50
xmin, xmax = 0, 1
nx, ny = 10, 50
mesh = RectangleMesh(Point(xmin, ymin), Point(xmax, ymax), nx, ny)
coordinates = mesh.coordinates()
V = FunctionSpace(mesh, 'CG', 1)

# Initial condition
#u_initial = interpolate(Expression("-x[1]", degree=1), V)
u_initial = Function(V)
u_initial.vector()[:] = 0.0

# Boundary conditions
## Dirichlet
class TopBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], ymax)

class BottomBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], ymin)

boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
TopBoundary().mark(boundaries, 1)
BottomBoundary().mark(boundaries, 2)
bc_bottom = DirichletBC(V, Constant(0.0), boundaries, 2)

## Neumann
q = Constant(1E-4) # as a source in the variational equation

# Variational problem
dt_ = 1000.0
T = 10000.0

## To store the previous solution and solve the next
u_old = Function(V)  # function to hold the value from the previous iteration

## To export the results
u_export = Function(V)

## To solve the problem
u = Function(V)
v = TestFunction(V)

## Elevation
z_expression = Expression('x[1]',degree=1)
z = project(z_expression, V)

## Variational equation
F = (theta(u) - theta(u_old))/dt_ * v * dx + K(u) * dot(grad(u + z), grad(v)) * dx + q * v * ds(1)

# Solver
## Parameters
solver = NonlinearVariationalSolver(NonlinearVariationalProblem(F, u, bc_bottom))
prm = solver.parameters
prm['newton_solver']['linear_solver'] = 'gmres'
prm['newton_solver']['maximum_iterations'] = 500
prm['newton_solver']['absolute_tolerance'] = 1E-8
prm['newton_solver']['relative_tolerance'] = 1E-7
prm['newton_solver']['relaxation_parameter'] = 0.8
prm['newton_solver']['report'] = True

## Time-stepping and data extraction
u_all_timesteps = []
y_coords_unique = np.unique(mesh.coordinates()[:, 1])
y_indices = np.where(np.isclose(y_coords_unique, 0.5))[0]

## Solving loop
t = 0
while t <= T:
    solve(F == 0, u, bc_bottom)
    t += dt_

# Plot data
x_center = (coordinates[:, 0].max() + coordinates[:, 0].min()) / 2.0
y_positions = np.unique(coordinates[:, 1])

u_at_center = []
for t in range(int(T/dt_)):
    u_values_for_t = []
    for y in y_positions:
        values = np.array([0.0])
        u_all_timesteps[t].eval(values, np.array([x_center, y]))

u_at_center = np.array(u_at_center)

for t in range(u_at_center.shape[0]):
    plt.plot(u_at_center[t, :], y_positions)

I’ve played around with different solver parameters, and tried a Picard approach, without success. I suspect issues with the variational form. I would deeply appreciate any insights or suggestions you might have.

Thank you for your time and assistance, and best regards,


Your variational form has integral over the domain for the source term, but you’re using ds(1) (ds is not defined anywhere in your code). I don’t know if it’s your mistake in the variational form or script. Anyway, you should create an example where people can just copy paste your code and run. Right now, your code seems incomplete to me.

1 Like

Hi Kei,

I edited the code with the imports, which I indeed forgot.

import numpy as np
from fenics import *
import matplotlib.pyplot as plt

The ds object is an integral type from the Measure class in ufl.

Thanks for your time,


You’re right about ds, my bad! But my point is that your variational form is \int_{\Omega} q v dx while ds is a surface integral, so there’s mismatch between your equation and script.

EDIT: signs on the second and third terms on the left-hand side.

My variational formulation was wrong. I edited my question consequently, so that the new equation is

\int_{\Omega} \frac{\partial \theta(u)}{\partial t} v ~ d\Omega + \int_{\Omega} K(u) \nabla (u + z) \cdot \nabla v ~ d\Omega + \int_{\partial\Omega} v K(u) \nabla (u + z) \cdot n ~ ds = 0

The form in FEniCS contained sign errors, and should be

## Variational equation
F = (theta(u) - theta(u_old))/dt_ * v * dx + K(u) * dot(grad(u + z), grad(v)) * dx - q * v * ds(1)

The Neumann condition q doesn’t seem to affect the results.

Thanks for your time,


What is your exact problem right now? I don’t know about Richards equation, but your dt seems to be very large to me. Are you sure that your choice of dt does not introduce numerical instability?

Hi, this is a soil column with a water table at its base, submitted to an inflow at the top (water precipitations). I used long time steps to be able to observe the evolution of the simulation. I obtained this plot


… but I expect something like this


Thanks again for your time,