Geoscience Reference
In-Depth Information
Box 5.3 Solution of Partial Differential Equations: Some Basic Concepts
As noted in the main text, it is nearly always not possible to obtain analytical solutions to
the nonlinear differential equations describing hydrological flow processes for cases of real
interest in practical applications, such as rainfall-runoff modelling. The approximate numerical
solution of nonlinear differential equations is, in itself, a specialism in applied mathematics
and writing solution algorithms is something that is definitely best left to the specialist. A major
reason for this is that it is quite possible to produce solution algorithms that are inaccurate or
are not
consistent
with the original differential equation (i.e. the approximate solution does
not converge to the solution of the original equation as the space and time increments become
very small). It is also easy to produce solutions for nonlinear equations that are not
stable
.
A stable solution means that any small errors due to the approximate nature of the solution
will be damped out. An unstable solution means that those small errors become amplified, often
resulting in wild oscillations of the solution variable at adjacent solution nodes or successive
time steps. The aim of this box is to make the reader aware of some of the issues involved in
approximate numerical solutions and what to look out for in using a model based on one of
the numerical algorithms available.
The differential equations of interest to the hydrologist (such as the Richards equation of
Box 5.1) generally involve one or more space dimensions and time. An approximate solution
then requires a discretisation of the solution in both space and time to produce a grid of points
at which a solution will be sought for the dependent variable in the equation. Figure B5.3.1
shows some ways of subdividing a cross-section through a hillslope into a grid of solution
points or nodes in space and Figure B5.3.2 shows a regular discretisation into both time and
space increments,
t
and
x
, for a single spatial dimension.
Starting a numerical solution always requires a complete set of nodal values of the solu-
tion variable at time
t
0. The algorithm then aims to step the solution through time, using
time steps of length
t
, to obtain values of the variable of interest for all the nodes at each
time step. The most easily understood numerical approximation method is the
finite differ-
ence
method. The original SHE model (Section 5.2.2) uses a finite difference solution of the
surface and subsurface flow equations on a regular spatial grid. The differentials in the flow
equation are replaced directly by differences. For example, to calculate a spatial differen-
tial for a variable
at node
i
in the interior of a grid at time step
j
, one possible difference
approximation is:
=
∂
∂x
i
+1
,j
−
i
−1
,j
2
x
i,j
≈
(B5.3.1)
This form is called a “centred difference approximation”. It is a more accurate approximation
than the forward difference that is often used for the time differential:
∂
∂t
i,j
+1
−
i,j
i,j
≈
(B5.3.2)
t
The forward difference is convenient to use for any differential in time since, given the nodal
values of
at time step
j
, the only unknowns are the values at time step
j
1. A centered
approximation for a second-order spatial differential involving a spatially variable coefficient,
K
, may be written as
∂
∂x
+
K
∂
∂x
i
+1
,j
−
K
i,j
x
i,j
i,j
−
i
−1
,j
x
i,j
≈
−
x
K
i,j
i
+1
,j
−
2
i,j
−
i
−1
,j
≈
x
2
