Solving the Richards Equation

I am currently working on solving the Richards Equation using the Brooks Corey Model using an Algebraic approach. I was trying to use a finite, crank nicolson or newton approach but I was difficulty as this area is still new to me. If there is an approach you would recommend over any other I am open to suggestions. Right now, my goal is to solve for theta and to plot the different values of theta in depth across time until an equilibrium is reached ( meaning theta is barely changing). I would also be interested in solving it with matlab's pdepe but I am a little confused how.

Richards Equation:

Brooks-Corey Model

Characteristic:

ConductiviTy:

soilData = {
    'Sand',        [0.395, 1.76e-2, 12.1, 4.05]; 
    'Silt Loam',   [0.485, 7.20e-4, 78.6, 5.30]; 
    'Clay',        [0.482, 1.28e-4, 40.5, 11.4];
};

%% --- Simulation Parameters ---
L = 100;         % cm (soil depth)
dx = 1;          % spatial resolution (cm)
dt = 0.001;         % time step (s)
T = 10 * 60;       
theta_r = 0.05;  % residual water content

x = 0:dx:L; 
Nx = length(x);
time = 0:dt:T; 
Nt = length(time);

%% --- Loop Over Soil Types ---
for i = 1:size(soilData, 1)
    % Extract soil parameters
    name = soilData{i,1};
    [phi, Kh, psi_ae, b] = deal(soilData{i,2}(1), soilData{i,2}(2), ...
                                soilData{i,2}(3), soilData{i,2}(4));
    theta_s = phi;
    initial_theta0 = theta_s;

    % Initialize moisture profile
    
    %theta_0 = theta_r + (theta_s-theta_r) * exp(-x/50);
    %theta = theta_0;
    
    theta = theta_s * exp(-x/70);
    %theta(1) = theta_s;
    theta_all = zeros(Nt, Nx);
    theta_all(1, :) = theta;
   
    % --- Time Loop (Algebraic Form) ---
    for t = 2:Nt
 
        % Initialize Deirvatives
        dKh_th_dz  = zeros(1,Nx);
        dpsi_th_dz = zeros(1,Nx);
        dKhdpsi_dz = zeros(1,Nx);
        % Kh_th=ones(1,Nx);

        % Compute theta_star
        theta_star = (theta - theta_r) / (phi - theta_r);
        theta_star = clip(theta_star,1E-6,1);

        Kh_th  = Kh * theta_star.^(2*b+3);
        psi_th = psi_ae * theta_star.^(-b); 
        
        % -- Compute dKh/dz & dPsi/dz ---
        for j = 1:Nx-1
            
           dKh_th_dz(j) = (Kh_th(j+1)-Kh_th(j))/dx;
           dpsi_th_dz(j) = (psi_th(j+1)-psi_th(j))/dx;
           
        end

        % Calculate Kh*dPsi/dz
        Khdpsi = Kh_th .* dpsi_th_dz;

        % Calculate d(Kh*dPsi/dz)/dz
        for j = 1:Nx-1
            
           dKhdpsi_dz(j) = (Khdpsi(j+1)- Khdpsi(j))/dx;
           
        end

        % Calculate dtheta/dt
        dtheta_dt = dKh_th_dz + dKhdpsi_dz;

        % Update theta
        theta = theta + dtheta_dt * dt;
        theta = clip(theta, theta_r, theta_s);

        % Boundary Condition
        theta(1) = theta_s;

        theta_all(t,:) = theta;

    end

    % --- Plot 1: Moisture Profile Over Time ---
    figure;
    cmap = jet(10);
    time_indices = round(linspace(1, Nt, 10));
    hold on;
    for k = 1:length(time_indices)
        plot(theta_all(time_indices(k), :), x, 'Color', cmap(k,:), 'LineWidth', 1.5);
    end
    xlabel('\theta (Volumetric Water Content)');
    ylabel('Depth (cm)');
    title(sprintf('%s: Moisture Profile Over Time', name));
    legend(arrayfun(@(t) sprintf('t = %.2f s', time(t)), time_indices, 'UniformOutput', false), 'Location', 'East');
    set(gca, 'YDir', 'reverse');
    grid on;
end

[Walter Roberson]

Missing "end" statement.

[Isabella Marie]

i have updated the code to include the missing end statement

[Walter Roberson]

I altered the legend format slightly.

[Sam Chak]

I repositioned the legend so that it doesn't obstruct the rainbow-colored curves.

[Isabella Marie]

The code is running but my values of theta values are not being updated correctly and I am not sure if that has to do due an error in how the richards equation is being approached. This is why the plot is not looking like it should.

[Walter Roberson]

I am not familiar with the Richards equation, so I do not know how the plots are supposed to look !

Answer 1

[Rahul]

Hi Isabella,

I understand that you are working on simulating 1D water infiltration in MATLAB using the Richards Equation along with the Brooks-Corey soil model, and you are trying to ensure that the current implementation accurately captures the underlying physics.

Based on the provided snippet, it looks like you are using an explicit time-stepping approach where the change in water content 'θ' is calculated using the mixed form of the Richards Equation.

This approach is mathematically sound and the use of Brooks-Corey parameters for hydraulic conductivity 'K(θ)' and pressure head 'ψ(θ)' should work in ideal scenarios.

However, as the system becomes stiffer (e.g., sharp wetting fronts or highly nonlinear behavior), explicit schemes may suffer from numerical instability or require extremely small time steps. To address this, you can try integrating MATLAB’s built-in 'pdepe' solver. This solver is designed for 1D parabolic and elliptic PDEs, and it automatically handles time and space discretization, which can simplify your code significantly and offer improved numerical stability.

Here’s a basic structure for using 'pdepe' function with the Richards Equation:

function [c, f, s] = richards_eqn(x, t, u, dudx)
    theta = u;
    K = Ks .* theta.^((2 + 3*lambda)/lambda);
    psi = psi_b .* (theta.^(-1/lambda));
    dpsidx = gradient(psi, x);  % or just use dudx if modeling in terms of psi
    c = 1;                      % Coefficient of du/dt
    f = K .* dpsidx + K;        % Flux term
    s = 0;                      % Source/sink (e.g., root uptake)
end

To solve the PDE, you can define the initial condition and boundary conditions separately, and then call 'pdepe' function as shown below:

sol = pdepe(0, @richards_eqn, @initial_condition, @boundary_conditions, x, t);

This approach should eliminate the need for manually coding finite-difference loops and can make your simulation more scalable and maintainable, especially if you plan to run it for long time spans or across different soil types.

For more information regarding the usage of 'pdepe' function, you can refer to the following documentation link:

• Solve 1-D parabolic and elliptic PDEs using 'pdepe' function: https://www.mathworks.com/help/matlab/ref/pdepe.html