Source: incomplete course notes for MatLab solvers, Session 4.2, from the NNESMO website:
https://nnesmo2019.webspace.durham.ac.uk/practicals/session42/
Dr Simon A Mathias,
Department of Earth Sciences,
Durham University
A confined aquifer is bounded to the West by a lake and to the East by an impermeable fault zone.
The aquifer is homogenous and isotropic and characterised by a transmissivity of 800 and a storativity of 0.01.
The edge of the lake lies parallel to the fault zone and is separated by a distance of 1600 m.
Initially the lake level is 50 m AOD (metres above ordinance datum).
After a significant episode of snow melt in the mountains above, the water level in the lake is suddenly raised to 53 m AOD.
The governing equation for one dimensional groundwater flow in a (un)confined aquifer takes the form
scipy: 1.9.1 numpy: 1.22.3
where S [-] is the aquifer storativity (which relates to the bulk compressibility of the aquifer),
h [L] is the hydraulic head, t [T] is time, L [L] is distance and
Q is [m^2/s] is the volumetric flow rate per unit breadth of confined aquifer, found from
where T is the transmissivity of the aquifer (which relates to the bulk permeability of the aquifer).
Discharge per unit width of aquifer =
The relevant initial and boundary conditions can be written as follows:
= 0
= 0.1, 0.3, 0.5, 1, 3, 5 in days
T= 800
S= 0.01
L= 1600
= = 50
= 53
and
Volume of Confined Water for a change in head = where:
Volume = volume drained from a confined aquifer for a hydraulic head change, Δh, over an area, A (L³)
S = storativity (dimensionless)
A = area over which the head change occurs (L²)
Δh = change in head (L)
Discharge per unit width of aquifer - (Measured in M³ per Second)
Transmissivity - (Measured in M² per Second)
Use the solver 'odeint' to develop a one-dimensional finite difference model to estimate the hydraulic head within the aquifer.
Show the results as a plot of hydraulic head against distance, for the following times: 0.1, 0.3, 0.5, 1, 3, 5 and 7 days.
The above problem is an example of a partial differential equation (PDE). However, if we discretise the problem in space, the problem becomes a coupled set of ordinary differential equations (ODE) with respect to time. Here we will use odeint to solve the large set of coupled ODEs that results from the spatial discretisation. We will use finite difference for spatial discretisation. However, other methods such as finite elements and pseudospectral methods can also be used in a similar way.
Let us consider a set of discrete points on the x-axis: .
The corresponding set of hydraulic heads can be written as: h(x).
An approximation of flow per unit breadth, , can be found from
but note that these flow rates are defined at an alternative set of points: .
The resulting set of ODEs to be solved takes the form
Letting ∆x → 0 and ∆t → 0 in (5) leads to the one-dimensional groundwater flow equation,
R=recharge, for 0<x<a, t>0, (6) where we have introduced the hydraulic diffusivity k via .
Recall that the relationship between and is
Furthermore it can be shown that
Create a new script-file:
This older solver has some perks:
lsoda
from the FORTRAN library odepack.20
So far all we have done is written some comments explaining what some of the subsequent variables are and defined the model parameter values.
The next step is to determine the locations on the x-axis that we are going to solve for.
In the above variable list we have x and xB. It is planned that x and xB will be vectors containing and , respectively.
Another point to note is that and .
Add the following code to provide a discretisation for 20 equally spaced solution points
array([ 0., 8., 16., 24., 32., 40., 48., 56., 64., 72., 80., 88., 96., 104., 112., 120., 128., 136., 144., 152., 160.])
array([ 4., 12., 20., 28., 36., 44., 52., 60., 68., 76., 84., 92., 100., 108., 116., 124., 132., 140., 148., 156.])
Next we will add some code to define the times of interest and a vector containing the initial values for .
This is the solver which should be used for new code, according to scipy. We need to use solver/method 'LSODA', instead of the default 'RK45' for handling possible stiffness.
so the new solver can use the finite diff. method.
(20,)
Of course we need to write the ODE function as well.
This function calculates derivatives and contains also a function for the analytical solution.
array([52.76103298, 52.80304448, 52.82861781, 52.8462694 , 52.85939305, 52.86964332, 52.8779358 , 52.8848236 , 52.89066321, 52.89569609, 52.90009236, 52.90397583, 52.90743901, 52.91055264, 52.9133718 , 52.91594018, 52.91829287, 52.92045847, 52.92246053, 52.92431864])
0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | 16 | 17 | 18 | 19 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 52.761033 | 52.803044 | 52.828618 | 52.846269 | 52.859393 | 52.869643 | 52.877936 | 52.884824 | 52.890663 | 52.895696 | 52.900092 | 52.903976 | 52.907439 | 52.910553 | 52.913372 | 52.915940 | 52.918293 | 52.920458 | 52.922461 | 52.924319 |
1 | 52.292531 | 52.414431 | 52.489350 | 52.541334 | 52.580113 | 52.610472 | 52.635074 | 52.655535 | 52.672900 | 52.687879 | 52.700972 | 52.712545 | 52.722870 | 52.732157 | 52.740569 | 52.748235 | 52.755259 | 52.761726 | 52.767707 | 52.773258 |
2 | 51.851225 | 52.041288 | 52.160358 | 52.243853 | 52.306556 | 52.355871 | 52.395972 | 52.429410 | 52.457846 | 52.482416 | 52.503922 | 52.522953 | 52.539949 | 52.555249 | 52.569118 | 52.581766 | 52.593361 | 52.604042 | 52.613924 | 52.623101 |
3 | 51.451782 | 51.692572 | 51.847807 | 51.958393 | 52.042277 | 52.108711 | 52.163007 | 52.208459 | 52.247233 | 52.280818 | 52.310276 | 52.336388 | 52.359743 | 52.380795 | 52.399899 | 52.417337 | 52.433338 | 52.448090 | 52.461747 | 52.474438 |
4 | 51.104361 | 51.375391 | 51.556901 | 51.688950 | 51.790461 | 51.871603 | 51.938372 | 51.994556 | 52.042683 | 52.084508 | 52.121295 | 52.153979 | 52.183271 | 52.209718 | 52.233754 | 52.255723 | 52.275905 | 52.294530 | 52.311789 | 52.327840 |
Polynomial function, f(x): 2 4.13e-10 x Derivative, f(x)'= 8.261e-10 x When x=5 f(x)'= 6.608623263102572e-07
array([ 1. , 1.47368421, 1.94736842, 2.42105263, 2.89473684, 3.36842105, 3.84210526, 4.31578947, 4.78947368, 5.26315789, 5.73684211, 6.21052632, 6.68421053, 7.15789474, 7.63157895, 8.10526316, 8.57894737, 9.05263158, 9.52631579, 10. ])
array([ 40., 120., 200., 280., 360., 440., 520., 600., 680., 760., 840., 920., 1000., 1080., 1160., 1240., 1320., 1400., 1480.])
[ 40. 120. 200. 280. 360. 440. 520. 600. 680. 760. 840. 920. 1000. 1080. 1160. 1240. 1320. 1400. 1480. 1560.] [52.2554889 52.69707656 52.77634693 52.81457874 52.83816115 52.85455826 52.86680633 52.87640383 52.88418678 52.89066321 52.89616202 52.90090674 52.90505534 52.90872303 52.91199607 52.91494051 52.9176079 52.92003912 52.92226709 52.92431864] (20, 20)
dtype('float64')