# Case study 1: Diffusion of fluid pressure and seismicity below Mt. Hood

We will apply our new transient model to study the relation between fluid pressure and seismicity in the crust below an active volcano, Mt. Hood in Oregon, USA. We will follow a publication by Saar and Manga (2003). The central claim of this paper is that there appears to be a correlation between seasonal recharge and watertable changes and seismicity in the upper crust. Saar and Manga suggest that this may be due to the effect of fluid pressure on seismicity. Groundwater recharge increases the watertable, which in turn increases pore fluid pressure in the crust. High fluid pressure means a lower effective normal stress on fault planes, which makes faults more likely to fail and generate seismic activity. We will use our 1D model to model the seasonal change in fluid pressure in the crust, and will try to quantify the magnitude of changes at a depth of 4500 m, where seismic activity is thought to originate.

## Workflow

1. Experiment with the timestep size to get a model that is stable
2. Model the effect of an instantaneous increase in recharge and record how long it takes before it reaches the depth where seismic activity occurs (4500 m below the surface)
3. Adjust hydraulic conductivity and storativity to find the optimal values that correspond to the phase shift between the recharge events and the spike in seimsic activity (151 days)
4. Implement a periodic boundary condition to simulate seasonal recharge and calculate the pore-pressure change at depth due to seismic activity

## 1. Keeping the model stable

The solution that we are using here is a so-called explicit finite difference solution. This is, next to the steady-state solution, the easiest and most intuitive approach to solving partial differential equations and simulating physical processes. However, explicit solutions have one drawback, the solution can become numerically unstable at large timestep sizes ($\Delta t$). We can test this experimentally. Increase the value of the timestep `dt` and run your model again. Keep on increasing the timestep until you get weird results. Record the timestep at which this occurs. Bonus points for the most artistic figure of unstable model results.

The timestep value at which these solutions become numerically unstable is predictable and has been quantified by three mathematicians from Göttingen, Courant, Friedrichs and Lewy in 1928. The stability condition is therefore also called the CFL or Courant number, and follows:

\begin{equation}
    CFL = \dfrac{q \Delta t}{\Delta x}
\end{equation}

The numerical solution becomes unstable for values of CFL that exceed 1.

**Assignment 1: Run the model several times and adjust the timestep size until the model becomes unstable. Choose a timestep size that is still stable, but large enough for the model to finish relatively fast. What timestep did you choose? Make a figure of an unstable model result.**

## 2. Model the effects of an instantaneous recharge event

We will start by setting up our numerical model and simulating the effect of an instantaneous increase in hydraulic head at the surface. The seasonal recharge at Mt. Hood is approximately equal to a layer of water of 1.5 m. The volcanic rocks at the surface have a specific yield of approximately 0.15. This means that a recharge of 1.5 m will increase the watertable by 1.5/0.15 = 10 m.   

We will try to implement an instantaneous change in hydraulic head at the surface by adjusting the top boundary condition *after* running the steady-state solution first. Add a line where you change the value of the hydraulic head at the surface to 10 m. Note that the variable controlling hydraulic head is an array instead of a single number like in exercise 2. This is because we want to make a boundary condition that varies over time later on.

**Assignment 2 Model the effects of an instantaneous increase in hydraulic head. How long does it take for the change to reach depths of 4500 m, where seismic activity occurs?**

## 3. Adjusting storativity and hydraulic conductivity

Next we will try to make our model more realistic. Instead of a single number for storativity ``S`` calculate the storativity using the specific storage equation given in your lecture handout (lecture 7). Note that we do not directly model the change in watertable. All the nodes in our model are assumed to be located below the watertable, and therefore we can ignore specific yield and assume that storativity is equal to the specific storage.

**Assignment 3: Implement the full equation for storativity in the model code and replace the line that currently assigns a value to ``S``. Use a reasonable value for fluid density, fluid compressibility and the compressibility of the rock matrix. Describe which parameter values you use and which value of storativity this results in. (see equation 16 and the parameter values mentioned below this equation. Use 1000 kg/m3 as density).**

The observed peak in seismicity underneath Mt. Hood shows a delay of approximately 151 days when compared to streamflow peaks at this location. Streamflow can be regarded as an indicator of seasonal groundwater recharge by snowmelt. The offset of peak seismicity may be an indicator of the time it takes for the fluid pressure increase at the surface to affect fluid pressures at 4500 m depth. 

**Assignment 4: Calibrate your model by adjusting ``K`` until the observed delay in increase in hydraulic head at 4500 m compared to the surface equals the observed delay in peak seismicity of ~151 days. Make figures of the model results. Which parameter value for K did you use for the final calibrated model?**

Note that changing K and/or S may make the model unstable again in some cases. This can be fixed by decreasing the timestep size.

## 4. Adding periodic boundary conditions

So far we have modeled the effect of an instantaneous groundwater recharge event. A more realistic way to model the effects of seasonal recharge is to include a periodic change in hydraulic head at the top boundary. Saar and Manga (2003) approximate this by using a cosine function (see equation 8). Uncomment the line where a periodic boundary condition is assigned in this notebook (in the transient parameters box below) and complete the line. Add a new box below, copy-paste the line to calculate the boundary condition and inspect the variable by typing ``print(h0)`` in the next line. Or ``print(h0[:100])`` if you want to inspect the first 100 values for instance. Make sure that you are modelling an increase from 0 to 10 m by checking the min and max values of ``h0`` by typing ``h0.min()`` and ``h0.max()``.

**Assignment 5 Run the model with the new boundary condition and make a figure of the result. What is the seasonal change in hydraulic head at a depth of 4500 m? What is the corresponding change in pore pressure (Pa)? Is this difference high or low compared to the difference between hdyraostatic and lithostatic pressure at these depths? (note, see lecture handouts for lecture on compaction, or ask the dutch guy in front of the classroom)**

# The actual model code:

## import python modules


```python
import matplotlib

%matplotlib inline

import numpy as np
import matplotlib.pyplot as pl
```

## Function to calculate steady-state hydraulic head


```python
def run_steady_state_model(x, dx, K, u0, W, n_iter=20000):
    
    C = (W * dx**2) / K

    u_new = np.ones_like(x) * u0
    u_old = np.ones_like(x) * u0

    # iterative solution steady-state h

    for t in range(n_iter):
    # make sure you indent anything below a for loop

        # set bnd conditions
        # left bnd:
        u_new[0] = u0
        # right bnd:
        u_new[-1] = 0.5 * (C[-1] + u_old[-1] + u_old[-2])

        # middle nodes:
        u_new[1:-1] = 0.5 * (C[1:-1] + u_old[2:] + u_old[:-2])

        u_old = u_new.copy()

    return u_new
```

## Model parameters


```python
day = 24.0 * 60.0 * 60.0
year = 365 * 24 * 60 * 60.0

L = 10000.0
dx = 100.0

# top boundary condition for the initial steady-state model
h0 = 0

# hydraulic conductivity
K = 1e-6
```

## Set up arrays


```python
# calculate depth of each node
z = np.arange(0, L + dx, dx)

# source term, we keep this at 0 for this exercise:
W = np.zeros_like(z)
```

## Parameters for transient model runs


```python
# specify storage coefficient
# for unconfined groundwater flow this is the specific yield
# for confined this is the specific storage
porosity = 0.1
density = 1000.0
g = 9.81
S = 1e-4

# formulate storativity using the compressibilites of water and the rock matrix
# uncomment and complete the following lines
# compressibility_water = ....
# compressibility_rocks = ....
# fluid_density = 1000.0
# g = 9.81
# S = ... a function of porosity, compressibility_water, compressibility_rocks, g, density

# timestep size
dt = 5.0 * day

#total duration
duration = 2.0 * year

# calculate total number of timesteps
n_timesteps = int(duration / dt)
```

## run steady-state model

We will first run the steady state model. The steady state value of h will be used as an initial value for the transient model runs


```python
h_steady = run_steady_state_model(z, dx, K, h0, W)
```

## Set up parameters for the transient model


```python
# set up array that records timesteps:
time = np.arange(n_timesteps) * dt

# define array to store flux and the variable over time:
n_nodes = len(z)
h = np.zeros((n_timesteps, n_nodes))

# set the steady-state value of u as value for first timestep
h[0] = h_steady

# increase the hydraulic head in the upper node:
h0 = np.zeros(n_timesteps)

# uncomment the following line to change the hydraulic head at the top
# boundary for the transient model runs:
h0[:] = 10

# uncomment the next lines to add a periodic boundary condition
#period = 1 * year
#amplitude = 10.0
#h0 = ... add a variation of equation 8 in the Saar and Manga (2003) paper here, in python language. 
# use the variables time, period, amplitude
# use np.pi for the number pi, np.cos(x) for the cosine of a variable x (replace x with your own desired function or number), 
# and use np.sin() for sine. 
```

## run the transient model:


```python
for j in range(1, n_timesteps):
    
    # calculate the flux between nodes
    q = -K * (h[j-1, 1:] - h[j-1, :-1]) / dx

    # set specified variable value at the left-hand node:
    h[j, 0] = h0[j]

    # implement no-flow boundary condition at right-hand side:
    q_right = 0.0
    h[j, -1] = h[j-1, -1] + (dt/S) * (-(q_right - q[-1])/dx) + (dt / S) * W[-1]

    # update nodes in the middle:
    h[j, 1:-1] = h[j-1, 1:-1] + (dt/S)*(-(q[1:] - q[:-1])/dx) + (dt/S) * W[1:-1]

    # print results to screen each 100 timesteps
    if j / 1000 == j / 1000.0:
        print('time = ', ((j * dt) / year), ', min, max value of h = ', h[j].min(), h[j].max())

```

    time =  0.0136986301369863 , min, max value of h =  0.0 10.0
    time =  0.0273972602739726 , min, max value of h =  0.0 10.0
    time =  0.0410958904109589 , min, max value of h =  0.0 10.0
    time =  0.0547945205479452 , min, max value of h =  0.0 10.0
    time =  0.0684931506849315 , min, max value of h =  0.0 10.0
    time =  0.0821917808219178 , min, max value of h =  0.0 10.0
    time =  0.0958904109589041 , min, max value of h =  0.0 10.0
    time =  0.1095890410958904 , min, max value of h =  0.0 10.0
    time =  0.1232876712328767 , min, max value of h =  0.0 10.0
    time =  0.136986301369863 , min, max value of h =  0.0 10.0
    time =  0.1506849315068493 , min, max value of h =  0.0 10.0
    time =  0.1643835616438356 , min, max value of h =  0.0 10.0
    time =  0.1780821917808219 , min, max value of h =  0.0 10.0
    time =  0.1917808219178082 , min, max value of h =  0.0 10.0
    time =  0.2054794520547945 , min, max value of h =  0.0 10.0
    time =  0.2191780821917808 , min, max value of h =  0.0 10.0
    time =  0.2328767123287671 , min, max value of h =  0.0 10.0
    time =  0.2465753424657534 , min, max value of h =  0.0 10.0
    time =  0.2602739726027397 , min, max value of h =  0.0 10.0
    time =  0.273972602739726 , min, max value of h =  0.0 10.0
    time =  0.2876712328767123 , min, max value of h =  0.0 10.0
    time =  0.3013698630136986 , min, max value of h =  0.0 10.0
    time =  0.3150684931506849 , min, max value of h =  0.0 10.0
    time =  0.3287671232876712 , min, max value of h =  0.0 10.0
    time =  0.3424657534246575 , min, max value of h =  0.0 10.0
    time =  0.3561643835616438 , min, max value of h =  0.0 10.0
    time =  0.3698630136986301 , min, max value of h =  0.0 10.0
    time =  0.3835616438356164 , min, max value of h =  0.0 10.0
    time =  0.3972602739726027 , min, max value of h =  0.0 10.0
    time =  0.410958904109589 , min, max value of h =  0.0 10.0
    time =  0.4246575342465753 , min, max value of h =  0.0 10.0
    time =  0.4383561643835616 , min, max value of h =  0.0 10.0
    time =  0.4520547945205479 , min, max value of h =  0.0 10.0
    time =  0.4657534246575342 , min, max value of h =  0.0 10.0
    time =  0.4794520547945205 , min, max value of h =  0.0 10.0
    time =  0.4931506849315068 , min, max value of h =  0.0 10.0
    time =  0.5068493150684932 , min, max value of h =  0.0 10.0
    time =  0.5205479452054794 , min, max value of h =  0.0 10.0
    time =  0.5342465753424658 , min, max value of h =  0.0 10.0
    time =  0.547945205479452 , min, max value of h =  0.0 10.0
    time =  0.5616438356164384 , min, max value of h =  0.0 10.0
    time =  0.5753424657534246 , min, max value of h =  0.0 10.0
    time =  0.589041095890411 , min, max value of h =  0.0 10.0
    time =  0.6027397260273972 , min, max value of h =  0.0 10.0
    time =  0.6164383561643836 , min, max value of h =  0.0 10.0
    time =  0.6301369863013698 , min, max value of h =  0.0 10.0
    time =  0.6438356164383562 , min, max value of h =  0.0 10.0
    time =  0.6575342465753424 , min, max value of h =  0.0 10.0
    time =  0.6712328767123288 , min, max value of h =  0.0 10.0
    time =  0.684931506849315 , min, max value of h =  0.0 10.0
    time =  0.6986301369863014 , min, max value of h =  0.0 10.0
    time =  0.7123287671232876 , min, max value of h =  0.0 10.0
    time =  0.726027397260274 , min, max value of h =  0.0 10.0
    time =  0.7397260273972602 , min, max value of h =  0.0 10.0
    time =  0.7534246575342466 , min, max value of h =  0.0 10.0
    time =  0.7671232876712328 , min, max value of h =  0.0 10.0
    time =  0.7808219178082192 , min, max value of h =  0.0 10.0
    time =  0.7945205479452054 , min, max value of h =  0.0 10.0
    time =  0.8082191780821918 , min, max value of h =  0.0 10.0
    time =  0.821917808219178 , min, max value of h =  0.0 10.0
    time =  0.8356164383561644 , min, max value of h =  0.0 10.0
    time =  0.8493150684931506 , min, max value of h =  0.0 10.0
    time =  0.863013698630137 , min, max value of h =  0.0 10.0
    time =  0.8767123287671232 , min, max value of h =  0.0 10.0
    time =  0.8904109589041096 , min, max value of h =  0.0 10.0
    time =  0.9041095890410958 , min, max value of h =  0.0 10.0
    time =  0.9178082191780822 , min, max value of h =  0.0 10.0
    time =  0.9315068493150684 , min, max value of h =  0.0 10.0
    time =  0.9452054794520548 , min, max value of h =  0.0 10.0
    time =  0.958904109589041 , min, max value of h =  0.0 10.0
    time =  0.9726027397260274 , min, max value of h =  0.0 10.0
    time =  0.9863013698630136 , min, max value of h =  0.0 10.0
    time =  1.0 , min, max value of h =  0.0 10.0
    time =  1.0136986301369864 , min, max value of h =  0.0 10.0
    time =  1.0273972602739727 , min, max value of h =  0.0 10.0
    time =  1.0410958904109588 , min, max value of h =  0.0 10.0
    time =  1.0547945205479452 , min, max value of h =  0.0 10.0
    time =  1.0684931506849316 , min, max value of h =  0.0 10.0
    time =  1.082191780821918 , min, max value of h =  0.0 10.0
    time =  1.095890410958904 , min, max value of h =  0.0 10.0
    time =  1.1095890410958904 , min, max value of h =  0.0 10.0
    time =  1.1232876712328768 , min, max value of h =  0.0 10.0
    time =  1.1369863013698631 , min, max value of h =  0.0 10.0
    time =  1.1506849315068493 , min, max value of h =  0.0 10.0
    time =  1.1643835616438356 , min, max value of h =  0.0 10.0
    time =  1.178082191780822 , min, max value of h =  0.0 10.0
    time =  1.1917808219178083 , min, max value of h =  0.0 10.0
    time =  1.2054794520547945 , min, max value of h =  0.0 10.0
    time =  1.2191780821917808 , min, max value of h =  0.0 10.0
    time =  1.2328767123287672 , min, max value of h =  0.0 10.0
    time =  1.2465753424657535 , min, max value of h =  0.0 10.0
    time =  1.2602739726027397 , min, max value of h =  0.0 10.0
    time =  1.273972602739726 , min, max value of h =  0.0 10.0
    time =  1.2876712328767124 , min, max value of h =  0.0 10.0
    time =  1.3013698630136987 , min, max value of h =  0.0 10.0
    time =  1.3150684931506849 , min, max value of h =  0.0 10.0
    time =  1.3287671232876712 , min, max value of h =  0.0 10.0
    time =  1.3424657534246576 , min, max value of h =  0.0 10.0
    time =  1.356164383561644 , min, max value of h =  0.0 10.0
    time =  1.36986301369863 , min, max value of h =  0.0 10.0
    time =  1.3835616438356164 , min, max value of h =  3.5348800511017475e-36 10.0
    time =  1.3972602739726028 , min, max value of h =  5.313631692816152e-35 10.0
    time =  1.4109589041095891 , min, max value of h =  4.702562351399871e-34 10.0
    time =  1.4246575342465753 , min, max value of h =  3.0863830579969544e-33 10.0
    time =  1.4383561643835616 , min, max value of h =  1.6559108715242797e-32 10.0
    time =  1.452054794520548 , min, max value of h =  7.640322532973994e-32 10.0
    time =  1.4657534246575343 , min, max value of h =  3.1291528655245264e-31 10.0
    time =  1.4794520547945205 , min, max value of h =  1.162072933417117e-30 10.0
    time =  1.4931506849315068 , min, max value of h =  3.973726770181997e-30 10.0
    time =  1.5068493150684932 , min, max value of h =  1.2656787237550293e-29 10.0
    time =  1.5205479452054795 , min, max value of h =  3.788787813944579e-29 10.0
    time =  1.5342465753424657 , min, max value of h =  1.0735736265308414e-28 10.0
    time =  1.547945205479452 , min, max value of h =  2.896332787766082e-28 10.0
    time =  1.5616438356164384 , min, max value of h =  7.475654216806541e-28 10.0
    time =  1.5753424657534247 , min, max value of h =  1.8535377639074024e-27 10.0
    time =  1.5890410958904109 , min, max value of h =  4.4300831529537575e-27 10.0
    time =  1.6027397260273972 , min, max value of h =  1.0237148300631569e-26 10.0
    time =  1.6164383561643836 , min, max value of h =  2.2931490763115765e-26 10.0
    time =  1.63013698630137 , min, max value of h =  4.990710683381966e-26 10.0
    time =  1.643835616438356 , min, max value of h =  1.057412582053079e-25 10.0
    time =  1.6575342465753424 , min, max value of h =  2.185026456811923e-25 10.0
    time =  1.6712328767123288 , min, max value of h =  4.41057306624575e-25 10.0
    time =  1.6849315068493151 , min, max value of h =  8.709337287331064e-25 10.0
    time =  1.6986301369863013 , min, max value of h =  1.684581880519881e-24 10.0
    time =  1.7123287671232876 , min, max value of h =  3.1954260163651656e-24 10.0
    time =  1.726027397260274 , min, max value of h =  5.9506167157784476e-24 10.0
    time =  1.7397260273972603 , min, max value of h =  1.0889766431303269e-23 10.0
    time =  1.7534246575342465 , min, max value of h =  1.9601528822971494e-23 10.0
    time =  1.7671232876712328 , min, max value of h =  3.4732511286285386e-23 10.0
    time =  1.7808219178082192 , min, max value of h =  6.063026304830206e-23 10.0
    time =  1.7945205479452055 , min, max value of h =  1.0434147870082653e-22 10.0
    time =  1.8082191780821917 , min, max value of h =  1.7714286877602342e-22 10.0
    time =  1.821917808219178 , min, max value of h =  2.9686177050970946e-22 10.0
    time =  1.8356164383561644 , min, max value of h =  4.913549945320053e-22 10.0
    time =  1.8493150684931507 , min, max value of h =  8.036687946572415e-22 10.0
    time =  1.8630136986301369 , min, max value of h =  1.299614763804055e-21 10.0
    time =  1.8767123287671232 , min, max value of h =  2.0787847067469636e-21 10.0
    time =  1.8904109589041096 , min, max value of h =  3.2904162194893586e-21 10.0
    time =  1.904109589041096 , min, max value of h =  5.156040450766675e-21 10.0
    time =  1.917808219178082 , min, max value of h =  8.00153297845471e-21 10.0
    time =  1.9315068493150684 , min, max value of h =  1.2302096292176574e-20 10.0
    time =  1.9452054794520548 , min, max value of h =  1.8744888194840886e-20 10.0
    time =  1.9589041095890412 , min, max value of h =  2.8315548806160036e-20 10.0
    time =  1.9726027397260273 , min, max value of h =  4.2416860891387943e-20 10.0
    time =  1.9863013698630136 , min, max value of h =  6.303031500820405e-20 10.0


## Some figures:

### A figure of h vs depth:


```python
fig, panel = pl.subplots(1, 1)

#  show change in hydraulic head over time
for j in range(0, n_timesteps, int(n_timesteps/20)):
    if j == 0:
        label = 't=0'
    else:
        label = None
    panel.plot(h[j], z, label=label)

label = 't=%0.1f yr' % (duration/year)
panel.plot(h[-1], z, color='black', lw=1.0, label=label)
panel.legend(loc='upper left', fontsize='small')

panel.set_ylabel('Depth (m)')
panel.set_xlabel('Hydraulic head (m)')
panel.set_ylim(L, 0)

fig.savefig('simulated_h_vs_depth.png')
```

### A figure of h over time for a specific depth


```python
# add the depths that you would like to show in the figures here (m)
target_depths = np.array([4500])

# add colors for the different depths
target_colors = ['blue', 'orange']

# look up the node for the different target depths
target_nodes = (target_depths / dx).astype(int)

fig, panel = pl.subplots(1, 1)

for target_depth, target_node, color in zip(target_depths, target_nodes, target_colors):
    label = '-%0.0f m' % target_depth
    panel.plot(time/year, h[:, target_node], color=color, label=label)

panel.legend()
panel.set_xlabel('Time (years)')
panel.set_ylabel('Hydraulic head (m)')

fig.savefig('simulated_h_over_time.png')
```

# References

Courant, R, K Friedrichs, and H Lewy. 1928. “Über Die Partiellen Differenzengleichungen Der Mathematischen Physik.” Mathematische Annalen 100 (1): 32–74.

Saar, M. O. & Manga, M. Seismicity induced by seasonal groundwater recharge at Mt. Hood, Oregon. Earth Planet. Sci. Lett. 214, 605–618 (2003).

Note, you can find these publications using google scholar: https://scholar.google.com




```python

```

ba6579817e6bde161a293f35d135a67b1f146fab
ipynb
Jupyter Notebook
exercises/exercise_3_transient_flow/.ipynb_checkpoints/exercise_3a_pore_pressure_diffusion_and_seismicity-checkpoint.ipynb
ElcoLuijendijk/fluids_in_the_crust
c2cadb0a91e9f9ed62094ac5e796168fef0d1a3e
CC-BY-4.0
2021-01-12T19:08:16.000Z
2021-01-13T14:27:42.000Z
exercises/exercise_3_transient_flow/.ipynb_checkpoints/exercise_3a_pore_pressure_diffusion_and_seismicity-checkpoint.ipynb
ElcoLuijendijk/fluids_in_the_crust
c2cadb0a91e9f9ed62094ac5e796168fef0d1a3e
CC-BY-4.0
exercises/exercise_3_transient_flow/.ipynb_checkpoints/exercise_3a_pore_pressure_diffusion_and_seismicity-checkpoint.ipynb
ElcoLuijendijk/fluids_in_the_crust
c2cadb0a91e9f9ed62094ac5e796168fef0d1a3e
CC-BY-4.0
Qwen/Qwen-72B
1. YES
2. YES