Environmental Engineering Reference
In-Depth Information
if
h0 > H
phi0 = -phi(iref,jref) + K*H*h0 - 0.5*K*H*H;
else
phi0 = -phi(iref,jref) + 0.5*K*h0*h0;
% reference potential
end
hc = 0.5*H+(1/K/H)*(phi+phi0);
% head confined
hu = sqrt ((2/K)*(phi+phi0));
% head unconfined
phicrit = phi0 + 0.5*K*H*H;
% transition confined / unconfined
confined = (phi>=phicrit);
% confined / unconfined indicator
h = confined.*hc+~confined.*hu;
% head
The reference value
phi0
is computed in the first five lines. If the specified
reference head
h0
exceeds the aquifer depth
H
, the aquifer is confined, otherwise it
is unconfined. For those situations the reference potential is computed using
different formulae, which are obtained by solving (
14.13
) for
'
0
.
Then the heads
hc
and
hu
for the confined and the unconfined situation are
both calculated.
phicrit
is the critical potential value at which the aquifer
switches between confined to unconfined state, according to formula (
14.14
).
The
confined
array contains a 1 at every entry corresponding to a mesh node,
where the aquifer is confined, and a 0 at every entry that corresponds with
a mesh node where the aquifer is unconfined. In the final command the real
2D head array is computed. Each entry related to a location in the confined part
takes the
hc
value, and the
hu
value if that is related to a location in
the unconfined part of the aquifer. '
~
' is the negation operator, switching an
1-entry to 0 and vice versa.
It is useful to indicate to the user the state of the aquifer within the model-region.
In the demonstration example below, messages concerning the state of the aquifer
are displayed in the command window. The
display
command produces a mes-
sage in the command window, showing the text string specified in the brackets.
if
all(all(confined))
display ('aquifer confined');
else
if
all(all(~confined))
display ('aquifer unconfined');
else
display ('aquifer partially confined and unconfined');
end
end
if
any(any(h<0))
display ('aquifer falls partially dry');
h = max(0, h);
end
For a confined aquifer the
confined
array contains only 1s. Here we use the
MATLAB
®
all
command to question if all entries are 1. The command has to be
used twice, as it operates on columns first, i.e. for a 2D array the command
all(A)
produces a line vector with the same number of columns as
A
. The column has
a 1-entry for each column containing non-zero elements only. The second call
gathers all columns. Some other MATLAB
commands operate on the same idea.
®
