<jupyter_start><jupyter_text>Adding Multiple WellsThis 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. **Authors**- _<NAME>, Artesia Water 2018_- _<NAME>, Artesia Water 2019_- _<NAME>, Artesia Water 2020_Import the requisite modules<jupyter_code>import numpy as np
import pandas as pd
import pastas as ps
import matplotlib.pyplot as plt

ps.show_versions()<jupyter_output>Python version: 3.7.4 (default, Aug  9 2019, 18:34:13) [MSC v.1915 64 bit (AMD64)]
Numpy version: 1.18.1
Scipy version: 1.4.1
Pandas version: 1.0.0
Pastas version: 0.15.0b<jupyter_text>Load data from a Menyanthes fileMenyanthes 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.<jupyter_code>fname = '../data/MenyanthesTest.men'
meny = ps.read.MenyData(fname)<jupyter_output><empty_output><jupyter_text>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 informatio later in our WellModel.<jupyter_code># 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))
df = pd.DataFrame(distances, index = extraction_names, columns=['Distance to Observation well'])
df<jupyter_output><empty_output><jupyter_text>Then plot the observations, together with the diferent stresses in the Menyanthes file.<jupyter_code># 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, val in enumerate(meny.IN.items()):
    name, data = val
    data["values"].plot(ax=axarr[i+1])
    axarr[i+1].set_title(name)
plt.tight_layout(pad=0)<jupyter_output><empty_output><jupyter_text>Create a model with a separate StressModel for each extractionFirst 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.<jupyter_code>ml = ps.Model(ps.TimeSeries(meny.H['Obsevation well']['values'], name="heads"))<jupyter_output>INFO: Cannot determine frequency of series heads<jupyter_text>Get the precipitation and evaporation timeseries and round the index to remove the hours from the timestamps.<jupyter_code>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"<jupyter_output><empty_output><jupyter_text>Create a recharge stressmodel and add to the model.<jupyter_code>sm = ps.StressModel2([IN, IN2], ps.Gamma, 'Recharge')
ml.add_stressmodel(sm)<jupyter_output>INFO: Inferred frequency from time series prec: freq=D 
INFO: Inferred frequency from time series evap: freq=D<jupyter_text>Get the extraction timeseries.<jupyter_code>stresses = []
for name in extraction_names:
    stress = ps.TimeSeries(meny.IN[name]['values'], name=name, settings='well')
    # the unit is allreay daily, so we need te set to_daily_unit to None
    stress.update_settings(to_daily_unit=None)
    stresses.append(stress)<jupyter_output>INFO: Cannot determine frequency of series Extraction 1
INFO: Time Series Extraction 1: values of stress were transformed to daily values (frequency not altered) with: divide
INFO: Cannot determine frequency of series Extraction 2
INFO: Time Series Extraction 2: values of stress were transformed to daily values (frequency not altered) with: divide
INFO: Cannot determine frequency of series Extraction 3
INFO: Time Series Extraction 3: values of stress were transformed to daily values (frequency not altered) with: divide<jupyter_text>Add each of the extractions as a separate StressModel.<jupyter_code>for stress in stresses:
    sm = ps.StressModel(stress, ps.Hantush, stress.name, up=False)
    ml.add_stressmodel(sm)<jupyter_output>WARNING: User-provided name 'Extraction 1' contains illegal character  
INFO: Time Series Extraction 1 was sampled down to freq D with method timestep_weighted_resample
WARNING: User-provided name 'Extraction 2' contains illegal character  
INFO: Time Series Extraction 2 was sampled down to freq D with method timestep_weighted_resample
WARNING: User-provided name 'Extraction 3' contains illegal character  
INFO: Time Series Extraction 3 was sampled down to freq D with method timestep_weighted_resample<jupyter_text>Solve the model.<jupyter_code>ml.solve()<jupyter_output>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: There are observations between the simulation timesteps. Linear interpolation between simulated values is used.<jupyter_text>Visualize the resultsPlot the decomposition to see the individual influence of each of the wells.In this graph we see that the influence of extraction 1 is very smooth.<jupyter_code>ml.plots.decomposition();<jupyter_output><empty_output><jupyter_text>We can calculate the gain of each extraction (quantified as the effect on the groudnwater level of an exttraction of 1000 m3/d).We can see that the gain of extraction 1 is largest, even though this extraction is furthest away from the observation well. The model has difficulty in separating the influence of the three extractions.<jupyter_code>for i in range(3):
    name = extraction_names[i]
    sm = ml.stressmodels[name]
    p = ml.get_parameters(name)
    gain = sm.rfunc.gain(p) * 1000000 / 365.25
    print("{0}: gain = {1:.2f} m / milion m^3/year".format(name, gain))
    df.at[name, 'gain StressModel'] = gain<jupyter_output>Extraction 1: gain = -0.70 m / milion m^3/year
Extraction 2: gain = -0.17 m / milion m^3/year
Extraction 3: gain = -0.13 m / milion m^3/year<jupyter_text>Create a model with a WellModelBecause we want to make the result more fysically feasable, we model the three extractions with a WellModel. This WellModel takes into acount the distances from the three extractions to the observation well, and assumes constant geohydrological properties.First we delete the existing StressModels with the well-data.<jupyter_code>for i in range(3):
    name = extraction_names[i]
    ml.del_stressmodel(name)<jupyter_output><empty_output><jupyter_text>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$, $\rho$ and $cS$, where parameter $\rho = \tfrac{r}{\lambda}$. This special version accepts 4 parameters: it interprets that fourth parameter as the distance $r$, and uses it to scale the $A$ and $\rho$ parameters accordingly. Create the WellModel and add to the model<jupyter_code>w = ps.WellModel(stresses, ps.HantushWellModel, "Wells", distances, settings="well")
ml.add_stressmodel(w)<jupyter_output>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<jupyter_text>Solve the model.As we can see, the fit with the measurements (EVP) is a litte lower now. The model is more fysically feasable however, because of the regularisation introduced by the WellModel.<jupyter_code>ml.solve()<jupyter_output>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
INFO: There are observations between the simulation timesteps. Linear interpolation between simulated values is used.<jupyter_text>Visualize the resultsPlot the decomposition to see the individual influence of each of the wells<jupyter_code>ml.plots.decomposition();<jupyter_output><empty_output><jupyter_text>Plot the stacked influence of each of the individual extraction wells in the results plot<jupyter_code>ml.plots.stacked_results(figsize=(10, 10));<jupyter_output><empty_output><jupyter_text>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 allways 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.<jupyter_code>wm = ml.stressmodels["Wells"]
for i in range(3):
    p = wm.get_parameters(model=ml, istress=i)
    gain = wm.rfunc.gain(p) * 1000000 / 365.25
    name = wm.stress[i].name
    print("{0}: gain = {1:.2f} m / milion m^3/year".format(name, gain))
    df.at[name, 'gain WellModel'] = gain<jupyter_output>Extraction 2: gain = -0.23 m / milion m^3/year
Extraction 3: gain = -0.17 m / milion m^3/year
Extraction 1: gain = -0.06 m / milion m^3/year<jupyter_text>Compare individual StressModels and WellModel Compare the gains that were calculated by the individual StressModels and the WellModel.<jupyter_code>df<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.07889490173739676, 'fraction_numerical': 0.016709389537643598, 'mean_word_length': 5.667088607594937, 'pattern_counts': {'":': 0, '<': 70, '<?xml version=': 0, '>': 70, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '19', 'dir': 'jupyter-structured-clean-dedup', 'id': '43720', 'license': 'MIT', 'n_tokens_mistral': 3342, 'n_tokens_neox': 3062, 'n_words': 1226, 'path': 'examples/notebooks/11_WellModel.ipynb', 'repo_name': 'RDWimmers/pastas', 'size': '10555'}