Modified Picard Iteration for the Mixed-Form Richards Equation

function [theta, K, cap] = Campbell (params, psi)

% -----------------------------
% Campbell (1974) relationships
% -----------------------------

% --- Soil parameters

theta_sat = params(1);    % Volumetric water content at saturation
psi_sat = params(2);      % Matric potential at saturation
b = params(3);            % Exponent
Ksat = params(4);         % Hydraulic conductivity at saturation

% --- Volumetric soil moisture (theta) for specified matric potential (psi)

if (psi <= psi_sat)
   theta = theta_sat * (psi / psi_sat)^(-1/b);
else
   theta = theta_sat;
end

% --- Hydraulic conductivity (K) for specified matric potential (psi)

if (psi <= psi_sat)
   K = Ksat * (theta / theta_sat)^(2*b+3);
else
   K = Ksat;
end

% --- Specific moisture capacity (cap) for specified matric potential (psi)

if (psi <= psi_sat)
   cap = -theta_sat / (b * psi_sat) * (psi / psi_sat)^(-1/b-1);
else
   cap = 0;
end

function [psi] = matric_potential (type, params, theta)

% --- Calculate psi for a given theta

switch type
   case 'van_Genuchten'

   theta_res = params(1);    % Residual water content
   theta_sat = params(2);    % Volumetric water content at saturation
   alpha = params(3);        % Inverse of the air entry potential
   n = params(4);            % Pore-size distribution index
   m = params(5);            % Exponent

   Se = (theta - theta_res) / (theta_sat - theta_res);
   psi = -((Se^(-1/m) - 1)^(1/n)) / alpha;

   case 'Campbell'

   theta_sat = params(1);    % Volumetric water content at saturation
   psi_sat = params(2);      % Matric potential at saturation
   b = params(3);            % Exponent

   psi = psi_sat * (theta / theta_sat)^-b;

end

function [soil] = picard (soil, params, dt, dpsi_tolerance, water_balance_error)

% -------------------------------------------------------------
% Use modified Picard iteration to solve the Richards equation
% -------------------------------------------------------------

% Input
% dt                   ! Time step (s)
% dpsi_tolerance       ! Convergence criterion for dpsi
% water_balance_error  ! Water balance error tolerance
% soil.nsoi            ! Number of soil layers
% soil.functions       ! van Genuchten or Campbell relationships
% soil.dz_plus_onehalf ! Thickness between between z(i) and z(i+1) (cm)
% soil.dz              ! Soil layer thickness (cm)
% soil.psi0            ! Soil surface matric potential boundary condition (cm)
%
% Input/output
% soil.theta           ! Volumetric soil moisture
% soil.psi             ! Matric potential (cm)
%
% Output
% soil.K               ! Hydraulic conductivity (cm H2O/s)
% soil.cap             ! Specific moisture capacity (/cm)
% soil.Q0              ! Infiltration flux (cm H2O/s)
% soil.QN              ! Drainage flux (cm H2O/s)
% soil.dtheta          ! Change in soil moisture (cm H2O)
% soil.err             ! Water balance error (cm H2O)

% --- Initialization

for i = 1:soil.nsoi

   % Save current soil moisture for time n
   theta0(i) = soil.theta(i);

   % Initialize delta_psi to a large value
   dpsi(i) = 1.0e36;

end

% --- Iteration loop

m = 0;
while (max(abs(dpsi)) > dpsi_tolerance) 

   % Increment iteration counter

   m = m + 1;

   % Stop if too many iterations

   if (m > 50)
      error ('Too many iterations')
   end

   % Hydraulic properties for current psi
   % theta - volumetric soil moisture
   % K -  hydraulic conductivity
   % cap - specific moisture capacity

   for i = 1:soil.nsoi
      switch soil.functions
         case 'van_Genuchten'
         [soil.theta(i), soil.K(i), soil.cap(i)] = van_Genuchten (params, soil.psi(i));
         case 'Campbell'
         [soil.theta(i), soil.K(i), soil.cap(i)] = Campbell (params, soil.psi(i));
       end
    end

   % Hydraulic conductivity at i+1/2 interface between layers i and i+1 is the arithmetic mean

   for i = 1:soil.nsoi-1
      K_plus_onehalf(i) = 0.5 * (soil.K(i) + soil.K(i+1));
   end

   % Hydraulic conductivity at i=1/2 between surface (i=0) and first layer i=1

   K_onehalf = soil.K(1);

   % dz at i=1/2 between surface (i=0) and first layer i=1

   dz_onehalf = 0.5 * soil.dz(1);

   % Terms for tridiagonal matrix

   i = 1;
   a(i) = 0;
   c(i) = -K_plus_onehalf(i) / soil.dz_plus_onehalf(i);
   b(i) = soil.cap(i) * soil.dz(i) / dt - a(i) - c(i);
   d(i) = K_onehalf / dz_onehalf * (soil.psi0 - soil.psi(i)) ... 
        - K_plus_onehalf(i) / soil.dz_plus_onehalf(i) * (soil.psi(i) - soil.psi(i+1)) ...
        + K_onehalf - K_plus_onehalf(i) - (soil.theta(i) - theta0(i)) * soil.dz(i) / dt;

   for i = 2:soil.nsoi-1
      a(i) = -K_plus_onehalf(i-1) / soil.dz_plus_onehalf(i-1);
      c(i) = -K_plus_onehalf(i) / soil.dz_plus_onehalf(i);
      b(i) = soil.cap(i) * soil.dz(i) / dt - a(i) - c(i);
      d(i) = K_plus_onehalf(i-1) / soil.dz_plus_onehalf(i-1) * (soil.psi(i-1) - soil.psi(i)) ...
           - K_plus_onehalf(i) / soil.dz_plus_onehalf(i) * (soil.psi(i) - soil.psi(i+1)) ...
           + K_plus_onehalf(i-1) - K_plus_onehalf(i) - (soil.theta(i) - theta0(i)) * soil.dz(i) / dt;
   end

   i = soil.nsoi;
   a(i) = -K_plus_onehalf(i-1) / soil.dz_plus_onehalf(i-1);
   c(i) = 0;
   b(i) = soil.cap(i) * soil.dz(i) / dt - a(i) - c(i);
   d(i) = K_plus_onehalf(i-1) / soil.dz_plus_onehalf(i-1) * (soil.psi(i-1) - soil.psi(i)) ...
        + K_plus_onehalf(i-1) - soil.K(i) - (soil.theta(i) - theta0(i)) * soil.dz(i) / dt;

   % Solve for the change in psi

   [dpsi] = tridiagonal_solver (a, b, c, d, soil.nsoi);

   % Update psi

   for i = 1:soil.nsoi
      soil.psi(i) = soil.psi(i) + dpsi(i);
   end

end

% --- Check water balance

soil.Q0 = -K_onehalf / dz_onehalf * (soil.psi0 - soil.psi(1)) - K_onehalf;
soil.QN = -soil.K(soil.nsoi);

soil.dtheta = 0;
for i = 1:soil.nsoi
   soil.dtheta = soil.dtheta + (soil.theta(i) - theta0(i)) * soil.dz(i);
end

soil.err = soil.dtheta - (soil.QN - soil.Q0) * dt;
if (abs(soil.err) > water_balance_error)
   error ('Water conservation error')
end

function [u] = tridiagonal_solver (a, b, c, d, n)

% Solve for U given the set of equations R * U = D, where U is a vector
% of length N, D is a vector of length N, and R is an N x N tridiagonal
% matrix defined by the vectors A, B, C each of length N. A(1) and
% C(N) are undefined and are not referenced.
%
%     |B(1) C(1) ...  ...  ...                     |
%     |A(2) B(2) C(2) ...  ...                     |
% R = |     A(3) B(3) C(3) ...                     |
%     |                    ... A(N-1) B(N-1) C(N-1)|
%     |                    ... ...    A(N)   B(N)  |
%
% The system of equations is written as:
%
%    A_i * U_i-1 + B_i * U_i + C_i * U_i+1 = D_i
%
% for i = 1 to N. The solution is found by rewriting the
% equations so that:
%
%    U_i = F_i - E_i * U_i+1

% --- Forward sweep (1 -> N) to get E and F

e(1) = c(1) / b(1);

for i = 2: 1: n-1
   e(i) = c(i) / (b(i) - a(i) * e(i-1));
end

f(1) = d(1) / b(1);

for i = 2: 1: n
   f(i) = (d(i) - a(i) * f(i-1)) / (b(i) - a(i) * e(i-1));
end

% --- Backward substitution (N -> 1) to solve for U

u(n) = f(n);

for i = n-1: -1: 1
   u(i) = f(i) - e(i) * u(i+1);
end

function [theta, K, cap] = van_Genuchten (params, psi)

% ----------------------------------
% van Genuchten (1980) relationships
% ----------------------------------

% --- Soil parameters

theta_res = params(1);   % Residual water content
theta_sat = params(2);   % Volumetric water content at saturation
alpha = params(3);       % Inverse of the air entry potential
n = params(4);           % Pore-size distribution index
m = params(5);           % Exponent
Ksat = params(6);        % Hydraulic conductivity at saturation
ityp = params(7);        % Soil texture flag

% --- Effective saturation (Se) for specified matric potential (psi)

if (psi <= 0)
   Se = (1 + (alpha * abs(psi))^n)^-m;
else
   Se = 1;
end

% --- Volumetric soil moisture (theta) for specified matric potential (psi)

theta = theta_res + (theta_sat - theta_res) * Se;

% --- Hydraulic conductivity (K) for specified matric potential (psi)

if (Se <= 1)
   K = Ksat * sqrt(Se) * (1 - (1 - Se^(1/m))^m)^2;

   % Special case for Haverkamp et al. (1977) sand (ityp = 1) and Yolo light clay (ityp = 2)

   if (ityp == 1)
      K = Ksat * 1.175e6 / (1.175e6 + abs(psi)^4.74);
   end
   if (ityp == 2)
      K = Ksat * 124.6/ (124.6 + abs(psi)^1.77);
   end

else

   K = Ksat;

end

% --- Specific moisture capacity (cap) for specified matric potential (psi)

if (psi <= 0)
   num = alpha * m * n * (theta_sat - theta_res) * (alpha * abs(psi))^(n-1);
   den =  (1 + (alpha * abs(psi))^n)^(m+1);
   cap = num / den;
else
   cap = 0;
end