Flow Modeling - Environmental Modeling - page 228

The coefficient of the pressure derivative on the left side is evaluated based on
the retention curve formula. For the van Genuchten formulation ( 11.23 ) one
obtains:

n 1

@y

@c ¼ y s y r

ð Þnma ah ð Þ

(11.25)

n mþ 1

ð 1 þ ac j j Þ

The differential ( 11.24 ) is the so-called Richards 11 equation , which can also be
found in a slightly different notation:

@y
@c

@c

@t ¼ @ @z K ðÞ @c @z

1 (11.26)

The following M-file is an implementation of the Richards equation with the van
Genuchten formulation for the suction-saturation and the conductivity-saturation
relationships.

L = 200; % length [L]
s1 = 0.5; % infiltration velocity [L/T]
s2 = 0; % bottom suction head [L]
T = 4; % maximum time [T]
qr = 0.218; % residual water content
f = 0.52; % porosity
a = 0.0115; % van Genuchten parameter [1/L]
n = 2.03; % van Genuchten parameter
ks = 31.6; % saturated conductivity [L/T]

x = linspace(0,L,100);
t = linspace(0,T,25);

options = odeset('RelTol',1e-4,'AbsTol',1e-4,'NormControl','off',…
'InitialStep',1e-7)
u=pdepe(0,@unsatpde,@unsatic,@unsatbc,x,t,options,... s1,s2,qr,f,a,n,ks);
figure;
title('Richards Equation Numerical Solution, computed with 100 mesh
points');

subplot (1,3,1);
plot (x,u(1:length(t),:));
xlabel('Depth [L]');
ylabel('Pressure Head [L]');

subplot (1,3,2);
plot (x,u(1:length(t),:)-(x'*ones(1,length(t)))');
xlabel('Depth [L]');
ylabel('Hydraulic Head [L]');
for j=1:length(t)
for i=1:length(x)
[q(j,i),k(j,i),c(j,i)]=sedprop(u(j,i),qr,f,a,n,ks);
end

end

subplot (1,3,3);
plot (x,q(1:length(t),:)*100)
xlabel('Depth [L]');
ylabel('Water Content [%]');

11 Lorenzo Adolph Richards (1904-1993), US-American soil physicist.