Issues on WellModel Scaling by a Distance

[quanerpython]

Describe the proposed feature
The WellModel Scaling by a Distance of the multiple_wells is likely underestimated other factors such as transmissivity. I have four extraction wells and got very good results with the stress model R2=95% but the results from the WellModel are not good, R2 =67% because the distance dominated the contribution to the head change at the observation well and the only closest extraction well contributes the head change, other the contribution of other three extraction wells is zero.

Expected usage
Quantify the percentage of the contribution from recharge, and each extraction well pumpage.
#Code to show the expected the behavior
The code should consider a parameter that scaling by both distance and transmissivity that was derived from the aquifer tests.
Additional context
The input data units are really confusing. I would like a details how to prepare the input files and their units associated. I would also suggest outputs of the percentage of each stress factor. The comparison plot of multiple models for my case is not readable and very difficult to modify. I would like to suggest an improvement for the comparison plot of the the multiple_wells.

[dbrakenhoff]

Hi @quanerpython,

Thanks for posting this issue! I'm converting it to a discussion, because I feel these are more questions to discuss, rather than direct issues with Pastas. But if any improvements to Pastas arise from this discussion I will open new Issues for those.

WellModel

First of all, the ps.HantushWellModel response function used by ps.WellModel assumes there is constant transmissivity. In fact, it assumes the system can be represented by a single confined aquifer with constant transmissivity. The only thing that affects the response is the distance between the pumping well and the observation well.

In your case, where you want more flexibility in the response and there are variations in transmissivity, you've already done exactly what we would recommend. In that case you split the pumping wells into multiple StressModels, each with their own response function. If you use ps.Hantush it still uses the same underlying approximations as ps.WellModel (single confining aquifer), but in this case the response (which is governed by, among other parameters, the transmissivity) can vary independently for each of the wells.

You mention aquifer tests in your question, what results from those tests would you want to include in your Pastas models?

Quantifying contributions

As for your question with regards to quantifying the contribution, this is challenging. One method is to compute the variance of each contribution and compare those values. However, with pumping, the variance of the contribution of the wells is highly dependent on the time period you select for your analysis. If the pump was already on, the variance of the contribution might be relatively small while the total drawdown might be significant. If you want to look at the absolute effect of each contribution and express it as a single statistic, how do you deal with negative and positive contributions? Another option is to compute the improvement of the model fit as you add stresses to it. This gives you some statistic of the relative importance of each stress. See for example this paper by Collenteur et al. (2023) (section 4.4). A final, simple option, that doesn't require too much thinking, might be to just plot the contribution of each stress over time in a single plot.

Input units

You are free to use any input units you want in Pastas. The only exception is the non-linear RechargeModel (ps.rch.FlexModel) which requires the precipitation and evaporation to be given in mm/day. The only reason we sometimes modify them is because we prefer to interpret the estimated model parameters in a certain way. For example,
the gain is given in Length / (unit of the stress). So for a pumping well this could be m / (m³/day), but for large pumping wells that might be a really small number, so then we translate the input stress to e.g. Mm³/year, so that the gain parameter is larger and more easy to analyze directly. These notebooks on the documentation website contain some information on units:

• Resampling time series
• Reasons to modify units

Perhaps we can improve the documentation further by including a clear section on units in one of the introductory notebooks. I will create an issue for this.

Comparing models

With lots of stresses it does become complex to create a clear comparison plot. In some cases you are better off creating your own comparison plot. We have tried to create some customizable comparison plots using the ps.CompareModels class. This notebook on the documentation website contains some detailed examples how you can customize your comparison. The trick with lots of wells is to plot them all in the same figure. That can be done using the smdict keyword argument, which allows you to list stresses that you want to include in the same subplot.

Hope this helps somewhat!

[quanerpython]

Hello @dbrakenhoff ,
Thank you so much for your reply and for taking further actions! I now have a clear understanding of the units for the input files. However, I believe it is essential to have a consolidated instruction on the input file format and the units of input and simulated parameters in one place.

I was able to plot individual outputs to avoid overlapping issues of multiple parameters in a single chart.

Here are more details about my testing project. Although all extractions likely tap into a hydraulically connected regional aquifer, the transmissivity varies spatially. My case is more general compared to assuming a uniform T for a region in present Pasta. The spatial variability of T is common for the alluvial aquifer due to coarse units like materials formed in paleochannels. Therefore, I believe it would be better to implement a well-specific T and consider the distance from the observation well.

////
The study area is situated in a desert region with minimal local recharge. There are five extraction wells, each located in separate sub-basins within a larger super-basin. These wells draw water from a regional carbonate aquifer that receives recharge from a distant northern region, with a travel time for the recharge estimated to exceed 100 years.

The overlying alluvial aquifer likely receives recharge from both artesian leakage from the regional carbonate aquifer and occasional local storm events. The extraction schedules are formulated based on a combined rate for all pumping wells within the basin, with each sub-basin's pumping centered on its respective extraction well.

Observational data and composite pumping schedules have been collected and are available for more than 20 years.

////

Specify coordinates observations

xo = 694305
yo = 4070248
The study area is located a desert region and local recharge is close to zero. The 5 extraction wells are
located five subbasins that are envloped within a superbasin.

The five exrcation wells are tapped the regional carbonate aquifer that receives its recharge from far
northern region. The traveling time from the recharge is likely over 100 years. The alluvial aquifer above
the regional carnobate aquifer is likley recharge by artisian leakage from the regional carbonate aquifer
plus some local recharge from local storm events.

The extraction schedules are prepared by a lumped rate of all pumping wells of the basin and assuming a
pumping center of Extraction Well for each basin.

The observations and the compoosite pumping schedules are availabe more than 20 years.

Specify coordinates and T extractions Extraction

WellModel_A			X		Y			T (ft2/d)
"Extraction_1": (689000, 4019000),    	80000
"Extraction_2": (688000, 4074000), 		250000 to 1000000      
"Extraction_3": (701000, 4068000),      310880   
"Extraction_4": (693000, 4037000),    	64000 to 530000
"Extraction_5": (697000, 4046000),     	80000

////////////////////////
Distance to observation well
Extraction_1 51618.619097
Extraction_2 7380.357850
Extraction_3 7242.653795
Extraction_4 33630.792851
Extraction_5 23781.857371

/////////////////////////////
Fit report Observation_well Fit Statistics

nfev 482 EVP 98.26
nobs 5397 R2 0.98
noise True RMSE 0.04
tmin 2006-05-03 00:00:00 AICc -39403.71
tmax 2023-04-11 00:00:00 BIC -39272.00
freq D Obj 1.81
warmup 3650 days 00:00:00 ___
solver LeastSquares Interp. No

Parameters (20 optimized)

                optimal     initial  vary

Recharge_A 0.049620 0.209252 True
Recharge_a 285.186705 10.000000 True
Recharge_f -1.526564 -1.000000 True
Extraction_1_A -0.000074 -0.001328 True
Extraction_1_a 3439.716578 100.000000 True
Extraction_1_b 0.005248 1.000000 True
Extraction_2_A -0.000052 -0.000157 True
Extraction_2_a 6772.820068 100.000000 True
Extraction_2_b 0.003974 1.000000 True
Extraction_3_A -0.000025 -0.000075 True
Extraction_3_a 7243.583773 100.000000 True
Extraction_3_b 0.006880 1.000000 True
Extraction_4_A -0.000291 -0.000326 True
Extraction_4_a 2599.352892 100.000000 True
Extraction_4_b 3.456402 1.000000 True
Extraction_5_A -0.000212 -0.001162 True
Extraction_5_a 3000.875861 100.000000 True
Extraction_5_b 0.384064 1.000000 True
constant_d 555.317005 553.247162 True
noise_alpha 5.059674 1.000000 True

Warnings! (5)

Response tmax for 'Extraction_1' > than calibration period.
Response tmax for 'Extraction_2' > than calibration period.
Response tmax for 'Extraction_3' > than calibration period.
Response tmax for 'Extraction_4' > than calibration period.
Response tmax for 'Extraction_5' > than calibration period.

///////////////////
Extraction_1: gain = -0.202 m / Mm^3/year
Extraction_2: gain = -0.141 m / Mm^3/year
Extraction_3: gain = -0.068 m / Mm^3/year
Extraction_4: gain = -0.796 m / Mm^3/year
Extraction_5: gain = -0.581 m / Mm^3/year

////////////////////
Fit report Observation_well Fit Statistics

nfev 160 EVP 91.96
nobs 5397 R2 0.92
noise True RMSE 0.10
tmin 2006-05-03 00:00:00 AICc -39006.37
tmax 2023-04-11 00:00:00 BIC -38947.06
freq D Obj 1.95
warmup 3650 days 00:00:00 ___
solver LeastSquares Interp. No

Parameters (9 optimized)

            optimal     initial  vary

Recharge_A 0.089091 0.209252 True
Recharge_n 1.756021 1.000000 True
Recharge_a 104.270116 10.000000 True
Recharge_f -1.829798 -1.000000 True
WellModel_A -2.537530 -0.000661 True
WellModel_a 869.264749 100.000000 True
WellModel_b -14.556620 -20.231614 True
constant_d 554.615014 553.247162 True
noise_alpha 24.317610 1.000000 True

Warnings! (1)

Parameter 'WellModel_b' on upper bound: -1.46e+01
///////////////////////////////