Table of Contents

One-dimensional transient groundwater flow

Scope

In this course, we will deal with transient groundwater flow in one-dimensional and radial situations (wells) for which analytic solutions are available. Analytic solutions are important because they allow insight in the behavior of the groundwater system, whereas numerical solutions do not; they only produce numbers. Analytic solutions are also important because they allow checking numerical models and checking numerical models is always necessary, not just because of possible errors in the model, but also because of possible errors in the input of the model. Analytical solutions also allow analysis of numerical models, which helps to understand their outcome. Finally, analytical solutions are powerful because they allow a rapid result with minimal input. They become even more powerful if combined with superposition and convolution.

Governing partial differential equation

We will always start our discussion with the governing differential equation at hand. Once we have it, we need to solve it. To be able to do that we need boundary conditions specifying fixed heads or fixed discharges along certain parts of the model boundaries. In the case of transient solutions we also need initial conditions specifying the head everywhere in the considered domain at time zero. Initial and boundary conditions are as important as the differential equation itself.

One-dimensional flow means a cross section with no flow components perpendicular to it.

We will treat analytical solutions for one layer only. Analytical solutions for more than one layer exist and have been extended to arbitrary numbers of layers in the 1980s by Kick Hemker and Kees Maas, see for instance Maas (1986) and Hemker and Maas (1987). These solutions require matrix computations, which were cumbersome at the time, but which may nowadays be readily computed in programs like Matlab and Maple. An implementation of the so-called eigen values and eigen vectors would be necessary to compute such solutions in a spreadsheet. This is possible, but beyond this course. Therefore, we limit ourselves to single-layer cases.

To derive the governing partial differential equation for the considered 1D situation, start with continuity. Consider a small section of the cross section. It has length $\Delta x$, the water balance over a short time $\Delta t$ may, continuity can be expressed as follows:

$$ \left( - \left( q_x + \frac {\partial q_x} {\partial x} \Delta x \right)\right)\,\Delta t + N \Delta t = S \left( \phi_{t + \Delta t} - \phi_t \right) $$

With $q$ m2/d specific discharge, $N$ recharge, $S$ the storage coeffiient of the aquifer and $\phi$ the head in the aquifer. $\Delta t$ is the considered time step and $\Delta x$ the considered spatial step.

Letting $\Delta x \rightarrow dx$ and $\Delta t \rightarrow dt$, we get

$$ - \frac {\partial q_x} {\partial x} + N = S \frac {\partial s} {\partial t} $$

with Darcy saying $q_x = -kD \partial \phi / \partial x$ we have

$$ kD \frac {\partial^2 \phi} {\partial x^2} + N = S \frac {\partial \phi} {\partial t} $$

We first considere the situation without recharge, so that $N=0$, giving

$$ kD \frac {\partial^2 \phi} {\partial x^2} = S \frac {\partial s} {\partial t} $$

Equation [eq:diffusion-equation-in-s] is known as the diffusion equation. It appears in many scientific fields like with diffusion, dispersion, heat conduction, sorption, consolidation and so on. Many researchers have derived solutions for it for specific boundary conditions. The coefficient $kD/S$ is called the diffusivity, often written as a thick $D$, like $\mathbb{D}$, which always has dimension [L2/T] whatever the scientific branch is. The diffusivity is the ratio of the ease of the flow (transmissivity) and the storage:

$$\mathbb{D}=\frac{kD}{S}$$

In the case of a phreatic (unconfined, water table) aquifer, the aquifer thickness is no longer constant. However, there are no transient solutions that take a time varying aquifer thickness into account. Linearization is then unavoidable, meaning that one has to choose a proper average aquifer thickness and be aware that the head change should remain small with respect to the saturated aquifer thickness.

Sinusoidal fluctuations in groundwater (like tides)

To compute flow and heads, we need to solve this partial differential equation with initial and boundary conditions pertaining to the situation we like to study.

This section deals exclusively with sinusoidal fluctuations in groundwater caused by a sinusoidal fluctuating head at the boundary at $x=0$. We deal with tidal fluctuations in groundwater first and then show temperature as a other application of the same basic partial differential equation.

A number of transient problems can be analyzed by assuming regular sinusoidal fluctuations of some boundary conditions. Generally the result will also behave like a sine in the same frequency. If we have the solution of a general harmonic fluctuation, we may solve many problems by superposition of any combination of frequencies that we need to adjust the fluctuations at the boundary to our liking. This way, hourly, daily, weekly and seasonal fluctuations may be readily combined. Examples of applications are tides in groundwater, penetration of temperature fluctuations, consolidation etc.

We will now consider a 1 dimensional aquifer perpendicular to a water body which is subject to sinus-like fluctuations of its water level. The situation is in the figure below.

Figure: Sinusoidal water level fluctuation in surface water causing tide in the groundwater system

The amplitude is indicated by the green arrow. The aquifer has constant transmissivity $kD$ and storage coefficient $S$. There is no recharge of the aquifer other than by interaction with the water body. The solution will be a sinusoidal fluctuation at any distance $x$, but it will be delayed and damped relative to the fluctuation in the surface water body. With this we try a solution to the pde that contains both damping with distance and delay.

The partial differential equation is

$$ \frac {\partial^2 \phi} {\partial x^2} = \frac S {kD} \frac {\partial \phi}{\partial t}$$

The presumed solution is

$$ s(x, t) = \phi(x, t) - \phi_0{x} = A \exp(-a x) \, \sin (\omega t - b x) $$

where $s(x, t)$ is the change of head due to the tidal fluctuations, which is independent of the mean head $\phi_0(x)$. Hence, $s(x, t)$ is superimposed on $\phi_0(x)$; it is irrelevant what $\phi_0(x)$ is.

$\omega$ is the angle velocity, such that $\omega T = 2 \pi$ with $T$ the cycle time (time of one complete wave).

The formula as written implies $s(0,0) = 0$ when $t=0$. One can choose any other angle at t=0, by writing the formula with

$$ \omega(t - t_0) = \omega\,t - \omega\,t_0 = \omega\,t -\beta_0 $$

or more conveniently, keep $\omega(t - t_0) - a x $ in the equation as it is.

Proof

To prove that this solution fulfills the partial differential equation, we have to fill it in and see if it is correct, and that it also fulfills any particular boundary and initial conditions.

So first derive its first and second derivatives with respect to $x$ and then its first derivative with respect to $t$. Because the amplitude is a constant, which will drop out, it can be left out of the derivation:

$$ \frac {\partial s} {\partial x} = -a e^{-a x} \sin(\omega t-b x)- b e^{-a x} \cos(\omega t- b x) $$$$ \frac {\partial^2 s} {\partial x^2} = a^2 e^{-a x} \sin(\omega t - b x) + a b e^{-a x} \cos(\omega t - b x) + a b e^{-a x} \cos(\omega t - b x) - b^2 e^{-a x} \sin(\omega t - b x)$$

$$ \frac S {kD} \frac {\partial s} {\partial t} = \frac S {kD}\omega e^{-a x} \cos(\omega t - b x)$$

By collecting the sines and cosines separately, we get

$$ a^2 -b^2 = 0 \,\,\,\, a = b $$

and, therefore,

$$ 2 a b = \frac S {kD}\omega $$

This yields

$$ a = \sqrt { \frac {\omega S} {2 kD} } $$

Note that this solution automatically fulfills the initial and boundary conditions.

Example, the intrusion of a continuous wave into the aquifer

Show the wave as it intrudes into the aquifer as it is agitated by the wave in the surface water at $x=0$. The values for this example are given in the code below.

It's nice to also indicate a point that surfs along with the wave

$$ \frac {\partial s}{\partial x} = - a e^{-a x}\left[\sin(\omega t - a x) + \cos(\omega t - ax)\right] = 0$$

$$ \sin(\omega t - a x) = - \cos(\omega t - a x) $$$$ \tan (\omega t - a x ) = -1 $$$$ \omega t - a x = - \frac \pi 4 + n \pi $$$$ x = \frac {\omega t - \frac \pi 4 + n \pi } a $$

Velocity of the wave

The graph shows the wave position for every hour and the red circles follow the top of the wave. The distance between two successive red points is the distance that the wave travels in 1 h. It is also the horizontal distance between a point at the same elevation at two successive waves. Verify that by inspecting the graph.

A formal derivation for the velocity of the top of the wave is given below, by computing the velocity of the point of the wave that has its first derivative equal to zero. It has the same answer: $v=\frac {dx} {dt} = \frac \omega \alpha$

$$\frac {\partial s} {\partial x} = -\alpha \exp(-\alpha x) \sin(\omega t-\alpha x) - \alpha \exp ( -\alpha x) \cos(\omega t- \alpha x) = 0 $$$$ \exp(-\alpha x) \left( \alpha \sin(\omega t - \alpha x) + \alpha \cos(\omega t - \alpha x) \right) = 0 $$$$ \exp(-\alpha x) \left( \tan(\omega t - \alpha x) + 1 \right) = 0 $$$$ tan(\omega t - \alpha x ) = -1 $$$$ \omega t - \alpha x = -\frac \pi 4 $$$$ \omega - \alpha \frac {dx} {dt} = 0 $$$$ \frac {dx} {dt} = \frac \omega \alpha $$

The semi-confined case (for illustration only)

In case we have a wave in a semi-confined aquifer with vertical leakage to or from an overlying layer with constant head, then no longer is $\alpha=\beta$. It is possible to derive the expressions, which then yields

$$ a = \frac 1 \lambda \sqrt { +\frac 1 2 + \frac 1 2 \sqrt{ 1 + \left( \omega S c \right)^2} } $$

and

$$ b = \frac 1 \lambda \sqrt{- \frac 1 2 + \frac 1 2 \sqrt{1 + \left( \omega S c \right)^2 } } $$

where $$ \lambda = \sqrt{ kD c} $$

It is easy to verity that we get when $ c \rightarrow \infty$, i.e. when the aquifer is confined.

What then happens is: The 1 under the deepest $\sqrt{ }$ diappears, and the popwer of 2 drops against the $\sqrt{}$. Then the first term $\frac 1 2$ under the top $\sqrt{}$ drops to give

$$ \alpha = \beta = \frac 1 \lambda \sqrt{\frac {\omega S c} 2 } = \sqrt {\frac {\omega S} {2 kD} } $$

which is the same as we derived for the confined situation.

Example

Now let us compare the semi-confined with a $\lambda = 600$ m, but further it is the same as the situation with the confined aquifer above.

Arbitrarily fluctuating head at x=0, superposition of waves

If we have an arbitrary fluctuating input signal, one that consists of several many simultaneous waves, each with its own amplitude, omega and phase shift, can we then compute its overal effect on the groundwater heads and flows?

Of course we can: this is just a matter of superposition, i.e. summing the effect of the individual waves.

Let us assume the fluctuation at $x=0$ can be described by 5 individual waves each with it's own amplitude, omega and phase shift. A phase shift is the phase at t=0. It is included in the argument of the $\sin(... - \theta)$ where $\theta = 2 \pi t_{delay}/T$ and $t_{delay}$ is the time at which the $\sin(\omega t -\theta) = 0$ at $x=0$.

Temperature fluctuations in the subsurface

For heat conduction the same diffusion equation applies if we replace head change by temperature. The only thing that changes is the Diffusivity. The ease of flow is now the heat conduction $\lambda$ [W/m] = [((E/T)/L$^2$)/(K/L)]=[E/(TKL)] and the storage is the heat capacity $\rho c$ [E/L$^3$/K]. The dimension is again [L$^2$/T]:

$$\mathbb{D} = \frac{\lambda}{\rho c} \mathtt{[E/(TKL) / E/(KL^3)]= [L^2/T]}$$

(notice that in the dimension E = energy, T = time, L = length, K = temperature (from Kelvin)).

Because both the heat conduction and the heat capacity have a contribution from both the water and the grains (solids) of the aquifer we can compute them as the combination of these contributions. With \epsilon for porosity we then have

$$\lambda = \epsilon\lambda_{w}+\left(1-\epsilon\right)\lambda_{s}$$$$\rho c = \epsilon\rho_{w}c_{w}+\left(1-\epsilon\right)\rho_{s}c_{s}$$

$\rho$ [M/L$^3$] is density and $c$ [E/(MK)], i.e. heat per kg solids per degree kelvin (= degree Celsius).

The heat capacity of saturated sandy soils is about $\lambda= 3$ W/m/K = 3 J/s/m/K. The specific heat capacity of water is $c_w$=4018 J/kg/K and that of sand grains $c_s \approx 800$ J/kg/K. With $\rho_s\approx 2650$ kg/m$^3$ and $\epsilon\approx 35\%$ we get $\rho c=2.85 \times 10^6$ J/m$^3$/K.

The diffusivity then becomes

$$\mathbb{D}=\frac{\lambda}{\rho c}=\frac{3}{2.85\times 10^6}=1.06 \times 10^{-3} \, \mathtt{m^2/s= 0.091 \, m^2/d}$$

From which we have

$$a=\sqrt{\frac{\omega}{2\mathbb{D}}}=\sqrt{\frac{\pi}{T\mathbb{D}}}$$

With these values it is possible to compute between what temperatures the temperature varies at different depth depending on the cycle time, i.e. due to daily, weekly, monthly or yearly temperature fluctuations at ground surface, the so-called temperature envelopes according to

$$T_{mean}-A\,\exp\left(-az\right)\le temp\le T_{mean}+A\,\exp\left(-az\right)$$

The monthly temperature in the Netherlands varies between 3.1C in January and 17.9 in July. The yearly amplitude is thus 7.4C with an average of 10.5C. Using this in the example, we can compute the temperature envelopes, that is, the lowest and highest temperatures lines between which the actual temperature will vary during one cycle time. Figure [fig:Temperature-envelopes-in-subsurface] shows the results, all computed for the same mean temperature and the same amplitude but for different cycle times as indicated in the legend. With the used data, the yearly temperature variation will barely reach 20 m below ground surface, a decade 50 m and 30 years 100 m. This implies that climate change may be measured using the change of the temperature at for instance 50 m below ground surface when measurements in the past are available, as is actually the case in the Netherlands.

Figure: Temperature envelopes in the subsurface due to temperature fluctuations at ground surface with mean 10$^o$C with different cycle times.

Illustration: An arbitrarily fluctuating river, solved by superposing its constituting waves using Fourier transforms

To make matters a bit more practical, let us assume that we have a measured time-varying river level for which we want to compute the dynamic effect on the groundwater, i.e. the time varying head and or flows.

So now not the sin or cos with its amplitudes, omegas and thetas are given, but the varying river level itself.

We could use our formula if we had this varying river level in the form of a number of individual waves, each of which has a fixed but arbitrary amplitude, frequency and phase.

So assuming that it is possible to split (decompose) a time-varying signal into a set of constituting waves, that, when added together give back that original varying river level, how can we do that?

It is possible to express any time varying signal into a large set of constituting waves. It is always done in telecommunications, audio, television and many other branches of science and technology. The tool for this is the Discrete Fourier Transform, or the Fast Fourier Transform (fft). It is readily available in python.

Let us do some exercises to show that it works.

First we construct a signal based on a set of known waves. Then we take this signal, try to find the underlying waves with the fast fourier transform, and then reconstruct the original signal.

Take a set of waves and compute the sum to be used as the signal

We have the (arbitrary) siginal, decompose it in its constituting waves using fft

As can be seen, the upper graph shows the frequencies that were discovered by the fft. The lower graph shows the original signal (blue, covered by the recovered signals. The first recovered signal used all frequencies, while the second recoverd signal only uses the frequencies with wave amplitude larger than 100. Both give essentially the same result.

The theory of the fft and the discrete fft can be found on Wikipedia. It is not necessary for the course, it is only used to demonstrate the techniquen and to show how a simple analytical groundwater solution can be utilized to simulate a complicated signal (or situation) by means of superposition.

Example with a measured river water level curve

Below is a "measured" varying water level of which we desired to compute the impact on the groundwater further away from the river.

The next step is to interpolate the curve and to sample it. Using scipy.interpolate.interp1d we generate a "function" that gives us interpolated values at any disired time.

Now let's see if we can break this signal into a set of waves for use by our function.

Which shows that we can reconstruct the the original signal accurataly with the waves that were discovered by fft.

This puts us in the position to construct the groundwater signal. This is done by superimposing the discovered waves.

Let's do an example.

Discrete Fourier transform (see Wikipedia)

Correlation between signal x and y is

$$ \sum_{i=0}^N x_i y_i $$

Correlation and DFT

$$ X_k =\sum_{i=0}^{N-1} x_n \exp \left( - 2 \pi i k \frac n N \right), \,\,\, k \in Z $$

This expresses the correlation between the signal $x(n)$ and the function $ \exp ( - 2 \pi i k \frac n N ) $

Back transformation:

$$ x_n = \frac 1 N \sum_{k=0}^{N-1} X_k \exp \left(2 \pi i k \frac n N \right), \,\,\,\, n \in Z $$

Power specttrum. This is in fact the absolute value of the fft outcome squared:

$$ P_k = Re \left(X_k^2\right) + Im \left(X_k^2 \right) = \left| X_k\right|^2$$

Amplitude is $ \left|X_x \right| $, the angle (phase) is $angle(X_x)$

$$ x_n = \sum_{k=0}^{N-1} \frac {X_k} N \left(\cos \left(\frac {2 \pi k n} N \right) - i \sin \left(\frac {2 \pi k n} N \right) \right) $$

Questions

1. Prove the correctness of the given solution yourself. As an extra exercise you could prove that equation [eq:QfluctuationStage] is correct by filling it into the partial differential equation for continuity $\frac{\partial Q}{\partial x}=-\frac{S}{kD}\frac{\partial s}{\partial t}$.

2. If the transmissivity if doubled, what is the effect on the drawdown?

3. When the \omega is doubled, what is the effect on the drawdown?

4. Explain how the distance to where the fluctuation of the sea or lake reaches in the aquifer depends on its frequency.

5. What is the depth effect of waves on the beach with one cycle per second, tides with one cycle per 12 hours, moon-tides with one cycle per two weeks in an aquifer with transmissivity is 500 m2/d and a storage coefficient of 0.001 and 0.1 respectively?

6. Let the solution to the diffusion equation for the confined aquifer be $s\left(x,t\right)=A\,\exp\left(-ax\right)\sin\left(\theta_{0}+\omega t-ax\right)$ and let $kD= 1000$ m2/d, $S=10e-3$, and the amplitude $A=2$ m, and $\theta_{0}$ is an arbitrary constant. Take time in days and show the head-change $s\left(x,t\right)$.

7. Include the discharge $Q\left(x,t\right)$

8. Compute and also show the envelope of the wave as a function of $x$.

9. How far inland can we still measure the tide if our device allows us to see a variation of 1 cm?

10. What is the velocity of the wave?

11. What is the delay at 1000 m from the shore (or show the delay graphically as function of $x$)?

12. Ad the case for a storage coefficient, $S_{y}=0.2$. And show the relation between the case with $S=0.001$ and $S_{y}=0.2$.

13. Create a complex input using of 4 sines, each with a different initial angle $\theta_{0}$, amplitude A and angular velocity $\omega$ and show the result.

Following question in a new computation:

1. The head in a lake above a clay bottom varies daily 30 cm. How deep does this fluctuation penetrate the underlying clay layer with conductivity of 10e-4 m/d and a specific storage coefficient of $S_{s}=0.0001$ /m?

2. What if the variation is weekly, monthly and seasonally only?

3. What is the penetration depth of a diurnal, seasonal and centennial temperature wave at ground surface given and porosity of $\epsilon=35%\%$ ?

4. If the sea is shallow and the clay layer is below the sea bottom, what will be the amplitude in the confined aquifer for the water compressibility $\beta_{w}=5\times10^{-5}$ and porous matrix compressibility $\beta_{s}=2e-5$ m2/N?

5. If the clay layer would be semi-pervious, what would this mean for the amplitude in the confined aquifer? Would it be greater, smaller compared with the case with a completely impervious layer?

Heat flow

1. Heat flow and groundwater flow show the same partial differential equation except that drawdown is replaced by temperature change. They both have one coefficient called diffusivity. What is the dimension of this diffusivity?

2. Diffusivity consists of a part expressing the ease of flow and a part that expresses the storage. What are these factors in the groundwater case and what are they in the heat flow case?

3. Show their dimensions and prove that the dimension of the diffusivity is the same in both cases.

4. Groundwater flow was ignored when we discussed heat flow. How would the temperature envelopes change due to upward groundwater seepage?

5. Also, how would they change due to downward seepage? How much would you estimate the effect on the yearly temperature envelopes if the recharge is 333 mm per year and porosity is 33\%?