Calibrate function - Accepting or rejecting results

[MattBrst]

I would like to use the following workflow:

1. I have a TTim model (ml1) with my initial parameters, and calculate heads (h_initial);
2. I make a Calibration run, based on some measurements, and recalculate heads (h_cal1)
3. I want to compare h_initial and h_cal1 to my measurements, and then accept or reject the parameters that were suggested by my calibration run.

So far I run into the following obstacles:

• The Calibrate run overwrites my initial parameters (or: assumes that I am happy with the result)
• I have not found an easy way to copy a model (something like: ml2 = ml1.copy()) which would help to preserve the original parameterset and experiment with another version of the model. Or is there a way to copy all parameters of the model and re-assign them somewhere?

Any help would be appreciated!

[mbakker7]

Thanks for the question, @MattBrst.
'deepcopy' is indeed what you want, but you to import it first from the copy module of the Python standard library.

import numpy as np
import ttim as ttm
from copy import deepcopy

# Data at Oude Korendijk taken from Kruseman and de Ridder (1990)
robs = 30 # distance to observation well, m
Q = 788 # discharge of well, m^3/d
time = np.array(
    [0.1 ,  0.25,  0.5 ,  0.7 ,  1.  ,  1.4 ,  1.9 ,  2.33,  2.8 ,
     3.36,  4.  ,  5.35,  6.8 ,  8.3 ,  8.7 , 10.  , 13.1 , 18.  ,
     27. ,  33. ,  41. ,  48. ,  59. ,  80. , 95.  , 139. , 181. , 
     245.,  300.,  360.,  480.,  600.,  728., 830.]) # in minutes
drawdown = np.array(
    [0.04 , 0.08 , 0.13 , 0.18 , 0.23 , 0.28 , 0.33 , 0.36 , 0.39 ,
     0.42 , 0.45 , 0.5  , 0.54 , 0.57 , 0.58 , 0.6  , 0.64 , 0.68 ,
     0.742, 0.753, 0.779, 0.793, 0.819, 0.855, 0.873, 0.915, 0.935,
     0.966, 0.99 , 1.007, 1.05 , 1.053, 1.072, 1.088]) # in meters
tobs = time / 24 / 60 # convert observation time to days
hobs = -drawdown # convert drawdown to heads

ml = ttm.ModelMaq(kaq=10, z=[-18, -25], Saq=1e-4, tmin=1e-5, tmax=1)
w = ttm.Well(ml, xw=0, yw=0, rw=0.2, tsandQ=[(0, Q)], layers=0)
ml.solve()

mlcal = deepcopy(ml) # copy of the model
cal = ttm.Calibrate(mlcal)  # Calibrate object
cal.set_parameter(name="kaq", initial=10, layers=0)  # Setting parameters
cal.set_parameter(name="Saq", initial=1e-4, layers=0)
cal.series(name="obs1", x=robs, y=0, t=tobs, h=hobs, layer=0)
cal.fit()

Now cal.rmse gives the root mean squared error of the calibrated model

cal.rmse()

give

0.031660181856832406

while computing the rmse of the original model still gives the rmse of the uncalibrated model

np.sqrt(np.mean((ml.head(robs, 0, tobs) - hobs) ** 2))

gives

1.5323521450225264

[MattBrst]

Could you elaborate the example a bit further? Let's take (an excerpt from) the WATEX example:

import numpy as np
import ttim

ml = ttim.ModelXsection(naq=1, tmin=1e-4, tmax=1e2)

z = [5, 0, -10]
kh = 10
c = 100
Saq = 1e-3
Sll = 1e-2
tsandhstar = [(0, 1)]

river = ttim.XsectionMaq(
    ml,
    -np.inf,  # river extends to infinitiy
    0,
    z=z,
    kaq=kh,
    Saq=Saq,
    Sll=Sll,
    c=c,
    topboundary="semi",
    tsandhstar=tsandhstar,
    name="river",
)

hinterland = ttim.XsectionMaq(
    ml,
    0,
    np.inf,  # hinterland extends to infinity
    kaq=kh,
    z=z,
    Saq=Saq,
    Sll=Sll,
    c=c,
    topboundary="semi",
    name="hinterland",
)

ml.solve()

ax = river.plot(labels=False, params=True, names=True)
hinterland.plot(ax=ax, params=True, names=True)
ml.elementlist[0].plot(ax=ax, hstar=10)  # hstar element river

from copy import deepcopy
mlcopy = deepcopy(ml)

So far so good. Now I want to assign values to the elements of the copied model. And see how it's changed. Something like:

river_copy.tsandhstar = [(0, -2), (10, 5)]
hinterland_copy.kaq = 10

ax = river_copy.plot(labels=False, params=True, names=True)
hinterland_copy.plot(ax=ax, params=True, names=True)
mlcopy.elementlist[0].plot(ax=ax, hstar=10)  # hstar element river

Would this be possible?

[mbakker7]

This is a very good question, but it should be a separate discussion, so that people can find the answer back. Can you start a new discussion with as title something like 'How to change time-varying boundary condition after model is created?'
Teaser: I have an answer.

[MattBrst]

Appreciated! New discussion: here. Would that also answer the second part of my question: How to change a hydrological parameter (kaq) of an existing model?

[MattBrst]

To answer my last question: Working further on the example of mbakker7/discussionmove#6, the following code works to change a parameter:

from copy import deepcopy
ml2 = deepcopy(ml)

# Check which items exist in your model:
for item in ml2.aq.inhomdict:
    print(item)

# Change parameter value kaq
ml2.aq.inhomdict['hinterland'].kaq = np.atleast_1d(10)

ml2.solve()

ax = ml2.aq.inhomdict['river'].plot(params=True, names=True)
ml2.aq.inhomdict['hinterland'].plot(ax=ax, params=True, names=True)
ax.set_ylim(ymax=10)