# Adding Multiple Wells¶

This notebook shows how a WellModel can be used to fit multiple wells with one response function. The influence of the individual wells is scaled by the distance to the observation point.

*Developed by R.C. Caljé, (Artesia Water 2020), D.A. Brakenhoff, (Artesia Water 2019), and R.A. Collenteur, (Artesia Water 2018)*

```
[1]:
```

```
import numpy as np
import pandas as pd
import pastas as ps
import matplotlib.pyplot as plt
ps.show_versions()
```

```
Python version: 3.7.8 | packaged by conda-forge | (default, Jul 31 2020, 02:37:09)
[Clang 10.0.1 ]
Numpy version: 1.18.5
Scipy version: 1.4.0
Pandas version: 1.1.2
Pastas version: 0.16.0b
Matplotlib version: 3.1.3
```

## Load data from a Menyanthes file¶

Menyanthes is timeseries analysis software used by many people in the Netherlands. In this example a Menyanthes-file with one observation-series is imported, and simulated. There are several stresses in the Menyanthes-file, among which are three groundwater extractions with a significant influence on groundwater head.

Import the Menyanthes-file with observations and stresses.

```
[2]:
```

```
fname = '../data/MenyanthesTest.men'
meny = ps.read.MenyData(fname)
```

Get the distances of the extractions to the observation well. Extraction 1 is about two times as far from the observation well as extraction 2 and 3. We will use this information later in our WellModel.

```
[3]:
```

```
# Get distances from metadata
xo = meny.H["Obsevation well"]['xcoord']
yo = meny.H["Obsevation well"]['ycoord']
distances = []
extraction_names = ['Extraction 1', 'Extraction 2', 'Extraction 3']
for extr in extraction_names:
xw = meny.IN[extr]["xcoord"]
yw = meny.IN[extr]["ycoord"]
distances.append(np.sqrt((xo-xw)**2 + (yo-yw)**2))
extraction_names = [name.replace(" ", "_") for name in extraction_names] # replace spaces in names for Pastas
df = pd.DataFrame(distances, index=extraction_names, columns=['Distance to Observation well'])
df
```

```
[3]:
```

Distance to Observation well | |
---|---|

Extraction_1 | 5076.464352 |

Extraction_2 | 2281.964490 |

Extraction_3 | 2783.783397 |

Then plot the observations, together with the diferent stresses in the Menyanthes file.

```
[4]:
```

```
# plot some series
f1, axarr = plt.subplots(len(meny.IN)+1, sharex=True, figsize=(14,7))
oseries = meny.H['Obsevation well']["values"]
oseries.plot(ax=axarr[0], color='k')
axarr[0].set_title(meny.H['Obsevation well']["Name"])
for i, (name, data) in enumerate(meny.IN.items(), start=1):
data["values"].plot(ax=axarr[i])
axarr[i].set_title(name)
plt.tight_layout(pad=0)
```

## Create a model with a separate StressModel for each extraction¶

First we create a model with a separate StressModel for each groundwater extraction. First we create a model with the heads timeseries and add recharge as a stress.

```
[5]:
```

```
ml = ps.Model(ps.TimeSeries(meny.H['Obsevation well']['values'], name="heads"))
```

```
INFO: Cannot determine frequency of series heads: freq=None. Resample settings are ignored and timestep_weighted_resample is used.
```

Get the precipitation and evaporation timeseries and round the index to remove the hours from the timestamps.

```
[6]:
```

```
IN = meny.IN['Precipitation']['values']
IN.index = IN.index.round("D")
IN.name = "prec"
IN2 = meny.IN['Evaporation']['values']
IN2.index = IN2.index.round("D")
IN2.name = "evap"
```

Create a recharge stressmodel and add to the model.

```
[7]:
```

```
sm = ps.StressModel2([IN, IN2], ps.Gamma, 'Recharge')
ml.add_stressmodel(sm)
```

```
INFO: Inferred frequency for time series prec: freq=D
INFO: Inferred frequency for time series evap: freq=D
```

Get the extraction timeseries.

```
[8]:
```

```
stresses = []
for name in extraction_names:
stress = ps.TimeSeries(meny.IN[name.replace("_", " ")]['values'], name=name, settings='well')
stresses.append(stress)
```

```
INFO: Cannot determine frequency of series Extraction_1: freq=None. Resample settings are ignored and timestep_weighted_resample is used.
INFO: Cannot determine frequency of series Extraction_2: freq=None. Resample settings are ignored and timestep_weighted_resample is used.
INFO: Cannot determine frequency of series Extraction_3: freq=None. Resample settings are ignored and timestep_weighted_resample is used.
```

Add each of the extractions as a separate StressModel.

```
[9]:
```

```
for stress in stresses:
sm = ps.StressModel(stress, ps.Hantush, stress.name, up=False)
ml.add_stressmodel(sm)
```

```
INFO: Time Series Extraction_1 was sampled down to freq D with method timestep_weighted_resample.
INFO: Time Series Extraction_2 was sampled down to freq D with method timestep_weighted_resample.
INFO: Time Series Extraction_3 was sampled down to freq D with method timestep_weighted_resample.
```

Solve the model.

Note the use of `ps.LmfitSolve`

. This is because of an issue concerning optimization with small parameter values in `scipy.least_squares`

. This is something that may influence models containing a WellModel (which we will be creating later) and since we want to keep the models in this Notebook as similar as possible, we’re also using `ps.LmfitSolve`

here.

```
[10]:
```

```
ml.solve(solver=ps.LmfitSolve)
```

```
INFO: Time Series Extraction_1 was sampled down to freq D with method timestep_weighted_resample.
INFO: Time Series Extraction_2 was sampled down to freq D with method timestep_weighted_resample.
INFO: Time Series Extraction_3 was sampled down to freq D with method timestep_weighted_resample.
INFO: Time Series Extraction_3 was extended to 1950-05-01 00:00:00 by adding 0.0 values.
INFO: There are observations between the simulation timesteps. Linear interpolation between simulated values is used.
```

```
Fit report heads Fit Statistics
========================================================
nfev 645 EVP 93.69
nobs 2844 R2 0.93
noise True RMSE 0.23
tmin 1960-04-28 12:00:00 AIC 19.97
tmax 2015-06-29 09:00:00 BIC 109.26
freq D Obj 22.83
warmup 3650 days 00:00:00 ___
solver LmfitSolve Interpolated Yes
Parameters (15 were optimized)
========================================================
optimal stderr initial vary
Recharge_A 1686.680335 ±19.54% 210.498526 True
Recharge_n 1.050298 ±3.71% 1.000000 True
Recharge_a 908.878934 ±27.43% 10.000000 True
Recharge_f -1.488694 ±14.46% -1.000000 True
Extraction_1_A -0.000872 ±3598.74% -0.000178 True
Extraction_1_a 9056.164653 ±958.60% 100.000000 True
Extraction_1_b 3.117841 ±2053.76% 1.000000 True
Extraction_2_A -0.000059 ±35.46% -0.000086 True
Extraction_2_a 2103.789349 ±203.32% 100.000000 True
Extraction_2_b 0.011679 ±236.92% 1.000000 True
Extraction_3_A -0.000021 ±39.92% -0.000170 True
Extraction_3_a 1069.202458 ±186.69% 100.000000 True
Extraction_3_b 0.003923 ±246.13% 1.000000 True
constant_d 11.141739 ±4.98% 8.557530 True
noise_alpha 49.583766 ±8.48% 1.000000 True
Parameter correlations |rho| > 0.5
========================================================
Recharge_A Recharge_a 0.86
Extraction_2_A -0.59
Recharge_n Recharge_a -0.68
Recharge_a Extraction_2_A -0.55
Recharge_f constant_d -0.99
Extraction_1_A Extraction_1_a 0.97
Extraction_1_b -1.00
Extraction_2_a -0.62
Extraction_2_b 0.60
Extraction_1_a Extraction_1_b -0.98
Extraction_2_a -0.60
Extraction_2_b 0.57
Extraction_1_b Extraction_2_a 0.67
Extraction_2_b -0.65
Extraction_2_A Extraction_2_a 0.84
Extraction_2_b -0.89
Extraction_2_a Extraction_2_b -0.99
Extraction_3_A Extraction_3_a 0.67
Extraction_3_b -0.84
Extraction_3_a Extraction_3_b -0.94
```

### Visualize the results¶

Plot the decomposition to see the individual influence of each of the wells.

```
[11]:
```

```
ml.plots.decomposition();
```

We can calculate the gain of each extraction (quantified as the effect on the groudnwater level of an extraction of ~1000 m3/d).

```
[12]:
```

```
for i in range(len(extraction_names)):
name = extraction_names[i]
sm = ml.stressmodels[name]
p = ml.get_parameters(name)
gain = sm.rfunc.gain(p) * 1e6 / 365.25
print("{0}: gain = {1:.2f} m / million m^3/year".format(name, gain))
df.at[name, 'gain StressModel'] = gain
```

```
Extraction_1: gain = -0.05 m / million m^3/year
Extraction_2: gain = -0.27 m / million m^3/year
Extraction_3: gain = -0.13 m / million m^3/year
```

## Create a model with a WellModel¶

We can reduce the number of parameters in the model by including the three extractions in a WellModel. This WellModel takes into acount the distances from the three extractions to the observation well, and assumes constant geohydrological properties. All of the extractions now share the same response function, scaled by the distance between the extraction well and the obervation well.

First we delete the existing StressModels with the well-data.

```
[13]:
```

```
for name in extraction_names:
ml.del_stressmodel(name)
```

We have all the information we need to create a WellModel: - timeseries for each of the extractions, these are passed as a list of stresses - distances from each extraction to the observation point, note that the order of these distances must correspond to the order of the stresses.

Note: the WellModel only works with a special version of the Hantush response function called `HantushWellModel`

. This is because the response function must support scaling by a distance \(r\). The HantushWellModel response function has been modified to support this. The Hantush response normally takes three parameters: the gain \(A\), \(a\) and \(b\). This special version accepts 4 parameters: it interprets that fourth parameter as the distance \(r\), and uses it to scale
the \(A\) and \(b\) parameters accordingly.

Create the WellModel and add to the model.

```
[14]:
```

```
w = ps.WellModel(stresses, ps.HantushWellModel, "Wells", distances, settings="well")
ml.add_stressmodel(w)
```

```
WARNING: It is recommended to use LmfitSolve as the solver when implementing WellModel. See https://github.com/pastas/pastas/issues/177.
INFO: Time Series Extraction_2 was sampled down to freq D with method timestep_weighted_resample.
INFO: Time Series Extraction_3 was sampled down to freq D with method timestep_weighted_resample.
INFO: Time Series Extraction_1 was sampled down to freq D with method timestep_weighted_resample.
```

Solve the model.

We are once again using `ps.LmfitSolve`

. The user is notified about the preference for this solver in a `WARNING`

when creating the WellModel (see above).

As we can see, the fit with the measurements (EVP) is the same as before.

```
[15]:
```

```
ml.solve(solver=ps.LmfitSolve)
```

```
INFO: Time Series Extraction_2 was sampled down to freq D with method timestep_weighted_resample.
INFO: Time Series Extraction_3 was sampled down to freq D with method timestep_weighted_resample.
INFO: Time Series Extraction_3 was extended to 1950-05-01 00:00:00 by adding 0.0 values.
INFO: Time Series Extraction_1 was sampled down to freq D with method timestep_weighted_resample.
INFO: There are observations between the simulation timesteps. Linear interpolation between simulated values is used.
```

```
Fit report heads Fit Statistics
======================================================
nfev 283 EVP 93.68
nobs 2844 R2 0.94
noise True RMSE 0.23
tmin 1960-04-28 12:00:00 AIC 8.06
tmax 2015-06-29 09:00:00 BIC 61.63
freq D Obj 23.02
warmup 3650 days 00:00:00 ___
solver LmfitSolve Interpolated Yes
Parameters (9 were optimized)
======================================================
optimal stderr initial vary
Recharge_A 1.296245e+03 ±17.40% 2.104985e+02 True
Recharge_n 1.002833e+00 ±3.65% 1.000000e+00 True
Recharge_a 8.534953e+02 ±28.45% 1.000000e+01 True
Recharge_f -2.000000e+00 ±13.14% -1.000000e+00 True
Wells_A -2.182820e-04 ±48.43% -8.609196e-05 True
Wells_a 6.959607e+02 ±32.20% 1.000000e+02 True
Wells_b 5.322245e-08 ±61.14% 8.749377e-08 True
constant_d 1.219999e+01 ±4.79% 8.557530e+00 True
noise_alpha 5.467031e+01 ±8.28% 1.000000e+00 True
Parameter correlations |rho| > 0.5
======================================================
Recharge_A Recharge_a 0.81
Wells_a -0.51
Recharge_n Recharge_a -0.72
Recharge_a Recharge_f -0.60
Wells_a -0.53
constant_d 0.77
Recharge_f constant_d -0.95
Wells_A Wells_a 0.93
Wells_b -1.00
Wells_a Wells_b -0.93
```

### Visualize the results¶

Plot the decomposition to see the individual influence of each of the wells

```
[16]:
```

```
ml.plots.decomposition();
```

Plot the stacked influence of each of the individual extraction wells in the results plot

```
[17]:
```

```
ml.plots.stacked_results(figsize=(10, 10));
```

Get parameters for each well (including the distance) and calculate the gain. The WellModel reorders the stresses from close to the observation well, to far from the observation well. We have take this into account during the post-processing.

The gain of extraction 1 is lower than the gain of extraction 2 and 3. This will always be the case in a WellModel when the distance from the observation well to extraction 1 is larger than the distance to extraction 2 and 3.

```
[18]:
```

```
wm = ml.stressmodels["Wells"]
for i in range(3):
p = wm.get_parameters(model=ml, istress=i)
gain = wm.rfunc.gain(p) * 1e6 / 365.25
name = wm.stress[i].name
print("{0}: gain = {1:.2f} m / million m^3/year".format(name, gain))
df.at[name, 'gain WellModel'] = gain
```

```
Extraction_2: gain = -0.23 m / million m^3/year
Extraction_3: gain = -0.17 m / million m^3/year
Extraction_1: gain = -0.04 m / million m^3/year
```

## Compare individual StressModels and WellModel¶

Compare the gains that were calculated by the individual StressModels and the WellModel.

```
[19]:
```

```
df
```

```
[19]:
```

Distance to Observation well | gain StressModel | gain WellModel | |
---|---|---|---|

Extraction_1 | 5076.464352 | -0.045164 | -0.044959 |

Extraction_2 | 2281.964490 | -0.273452 | -0.233410 |

Extraction_3 | 2783.783397 | -0.129116 | -0.169794 |