Solving a PDE problem considering the nonlinear diffusion coefficient of Richards equation

Question

I'm trying a very simple 1D simulation on Richards equation, simulating a soil column, with zero water pressure (Dirichlet boundary, u=0) at its base (left in 1D), and a rain (Neumann boundary, q=1E-5) at the top (right in 1D). In strong form,

The diffusion coefficient (k(u), or hydraulic conductivity) relies on water pressure (which should be negative all through the soil column due to capillarity). The variable z is the elevation from a reference point, here at the bottom (or left) of the mesh.

I coded the following script, but still can't figure out how to obtain results any close to the one I obtained using SEEP/W, a professional software solving Richards equation. I used a loop to recompute k(u) since it relies on u.

import fipy as fp
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# Simulation with professional software SEEP/W
seep = pd.read_csv('https://gist.githubusercontent.com/essicolo/e670981aa0171a449a554d7763157374/raw/8d06b07117f6f31a57b2156d07a2ba971a355e97/gistfile1.txt')

# Material parameters
alpha = 2.0 / 9.807  # [1/m]
n = 1.5
m = 1 - 1/n
l = 0.5
Ks = 1.0E-4  # [m/s]

# Top infiltration
rain = 5E-5

# Domain
nz = 200
Lz = 20.0
dz = Lz / nz
mesh = fp.Grid1D(nx=nz, dx=dz)

# Variables
z = fp.CellVariable(name="elevation", mesh=mesh, value=mesh.cellCenters[0])
u = fp.CellVariable(name="matrix pressure", mesh=mesh, value=-z.value)
u.setValue(u_min, where = z >= -u_min) # close to the real solution
u_old = fp.CellVariable(name="old matrix pressure", mesh=mesh, value=-u.value)
u_convergence = fp.CellVariable(name="convergence check", mesh=mesh, value=0.)
K = fp.CellVariable(mesh=mesh, value=Ks)

# Boundary conditions
u.constrain(0., mesh.facesLeft) # bottom
g = fp.FaceVariable(mesh=mesh, value=-rain)
u.faceGrad.constrain(g * mesh.faceNormals, mesh.facesRight) # top

# Solver
max_iters = 100
tolerance = 1e-10
solver = fp.LinearLUSolver()

convergence = []
iterations = []
for iter in range(max_iters):

    plt.plot(u.value, z.value, label=f'Iter {iter}')
    u_old.value = u.value

    K.setValue(Ks * (1 - (-alpha * u)**(n-1) * (1 + (-alpha * u)**n)**(-m)) * (1 + (-alpha * u)**n)**(-m*l), where=u < 0)
    K.setValue(Ks, where=u >= 0)
    
    eq = (fp.DiffusionTerm(coeff=K, var=u) + K == 0)
    eq.solve(var=u, solver=solver)

    u_convergence.value = np.abs(u.value - u_old.value)

    iterations.append(iter)
    convergence.append(np.mean(u_convergence.value))
    
    if np.all(u_convergence.value < tolerance):
        break


plt.plot(seep['u'], seep['z'], c='k', label='with SEEP/W')
plt.legend()

#plt.plot(iterations, convergence)

Thanks for your time and support,

Essi

Answer 1

fp.DiffusionTerm(coeff=K, var=u) + K == 0

isn't your equation.

Just as you have in 1D, in multiple dimensions:

where is the unit vector in the z direction.

I find that

eq = (fp.DiffusionTerm(coeff=K, var=u) + (K.faceValue * [[1]]).divergence == 0)

behaves more like your SEEP/W result, although it's hard to be sure without knowing what u_min is.