Calibration of loading efficiency values

[MattBrst]

Today I tried to calibrate on the loading efficiency with:

 cal.set_aquifer_parameter(
            "leffaq",
            layers=[1],
            initial=leffaq1_river,
            pmin=0,
            pmax=1,
            log_scale=True,
            inhoms=("river"),
            model="transient",
        )

But this results in: ValueError: Parameter 'leffaq' not recognised. Supported: ['kaq', 'c', 'Saq', 'Sll', 'kzoverkh']

So I guess we would need to add leffaq and leffll in this block:

timflow/timflow/calibrate.py

Lines 1381 to 1396 in 882bc67

	def _get_aquifer_parameter_array(model, aq, name: str) -> np.ndarray:
	"""Return reference to the appropriate parameter array in the aquifer object."""
	lookup = {"kaq": aq.kaq, "c": aq.c}
	if hasattr(aq, "Saq"):
	lookup["Saq"] = aq.Saq
	if hasattr(aq, "Sll"):
	lookup["Sll"] = aq.Sll
	if hasattr(aq, "kzoverkh"):
	lookup["kzoverkh"] = aq.kzoverkh
	assert aq.model == model, "Aquifer does not belong to the given model"
	for prefix, arr in lookup.items():
	if name.startswith(prefix):
	return arr
	raise ValueError(
	f"Parameter '{name}' not recognised. Supported: {list(lookup.keys())}"
	)

Or would you guys suggest something else?

[mbakker7]

Neither leffaq nor leffll are currently implemented as parameters that can be calibrated by specifying their string name.
The reason we have not done it (yet) is that these problems already often have correlated parameters that cannot be calibrated independently before the loading efficiency parameters were added. Luckily for the loading efficiency, the values don't vary over a wide range, and reasonable estimates can be made. For example leffaq=0.8 in sand layers or leffll=1 in a clay layer. My suggestion is to do a calibration w/o loading efficiency (leffaq=lefll=0) and then one by setting them to reasonable values. If the estimated aquifer parameters differ significantly, then the loading efficiency apparently plays an important role and you may investigate further.

Regarding estimating them, we can consider building them in so they can be estimated by name, like kaq or Saq. Right now, you can already estimate them by specifying them by reference, see this notebook: https://timflow.readthedocs.io/en/latest/transient/05benchmarks/synthetic_test_calibrate.html

[MattBrst]

Thanks for the explanation!

I got the example working with the suggested function:

cal.set_parameter_by_reference(name="leffaq", parameter=riv_tr.leffaq[:], initial=0.5, pmin=0.0, pmax=1.0, log_scale=True)

One minor thing worth noting: The function set_parameter_by_reference does not have the option to make a fit for two inhoms (let's say river and foreland) combined. To overcome this limitation, I would say it is more logical to add leffaq to set_parameter() than to add the coupled inhoms handling to set_parameter_by_reference().

[dbrakenhoff]

I don't think it would be very difficult to allow set_parameter_by_reference() to accept a list of arrays as parameter to deal with those cases.

I think adding a loop over the parameters when passed as a list or tuple, and adding each array separately as a calibration "target" for that parameter would allow it to work with an arbitrary collection of arrays across models or inhoms.

So this:

timflow/timflow/calibrate.py

Lines 499 to 510 in 50d7259

	assert isinstance(parameter, np.ndarray), "parameter must be a numpy array"
	if initial is None:
	initial = parameter[0]
	param = Parameter(
	name=name, initial=initial, pmin=pmin, pmax=pmax, log_scale=log_scale
	)
	param.add_target(parameter, slice(None))
	if log_scale:
	param.set(np.log10(np.abs(initial)))
	else:
	param.set(initial)
	self._parameters[name] = param

Becomes (including small fix check for non-zero initial when log_scale=True):

        param = Parameter(
            name=name, initial=initial, pmin=pmin, pmax=pmax, log_scale=log_scale
        )
        # allow collections of parameter arrays
        if isinstance(parameter, (list, tuple)):
            assert all(isinstance(ipar, np.ndarray) for ipar in parameter), (
                "parameter(s) must be a numpy array"
            )
            for ipar in parameter:
                param.add_target(ipar, slice(None))
            if initial is None:
                initial = parameter[0][0]  # first entry of first array
        # single parameter array
        else:
            assert isinstance(parameter, np.ndarray), "parameter must be a numpy array"
            param.add_target(parameter, slice(None))
            if initial is None:
                initial = parameter[0]
        if log_scale:
            if initial == 0:
                raise ValueError("initial must be nonzero when log_scale=True")
            param.set(np.log10(np.abs(initial)))
        else:
            param.set(initial)
        self._parameters[name] = param

Feel free to test this @MattBrst and submit a PR if the above works.

[MattBrst]

Thank you for the suggested code. I wrote another Q&A question #189 as I would say it would be more logical to distinguish between:

1. set_parameter() for all aquifer-related parameters (including the coupling of inhom params)
2. set_parameter_by_reference() for non-aquifer parameters (like well parameters)

If we proceed with your code above, I think we also need to check that pmin>0 for the logscale.

[mbakker7]

I will move this Issue to Discussions