Extract a Penman-Monteith submodel from HydPy-L(ARSIM)

[tyralla]

This is a specialised follow-up issue of #109, which serves as a notebook for our current and later efforts of optimising HydPy's Penman(-Monteith)-based evapo(transpi)ration methods.

[tyralla]

lland_v3 vs. lland_v4

We can target both application models in one step. This saves a few intermediate steps and should be no considerable complication, as both models rely on nearly identical calculation methods. The only difference that I am currently aware of is that lland_v4 suppresses Penman(-Monteith)-based interception evaporation from snow-covered tree canopies because the Knauf-based snow interception evaporation method steps in. lland_v3, which does not come with a snow interception module, calculates interception evaporation from trees always after Penman(-Monteith). "Trees" means land types LAUBW, MISCHW, and NADELW.

[tyralla]

actual albedo

lland_v3 and lland_v4 calculate net shortwave radiation based on time-variable albedo values. Without a snow cover, the actual albedo depends on the current month and land type. If there is a snow cover, the actual albedo depends on the snow's age. Unfortunately, this indicates a closer connection between the evaporation and the planned snow module than expected.

This is no problem for non-tree land types because snow here suppresses interception evaporation (EvI) and soil evapotranspiration (EvB) so that the actual albedo value becomes irrelevant for the AET submodel. However, for tree land types, lland_v3 and lland_v4 always calculate soil evapotranspiration (which seems not very consistent with the above observation on interception evaporation).

It seems the main model must be aware of the actual albedo. Before further thoughts, we should search for similar problematic connections.

[tyralla]

roughness length

LARSIM's default for calculating the roughness length is dividing the crop height by ten. However, lland_v3 and lland_v4 always use LARSIM's RAUHIGKEITSL. NACH QUAST UND BOEHM, 1997 option: $z_0 = 0.021 + 0.163 \cdot CropHeight$. LARSIM's documentation states that this option results in higher evapotranspiration estimates.

This is not a topic of modularisation. However, I mention it here because we should probably also support selecting different approaches for calculating the roughness length.

[tyralla]

soil surface resistance

For soils with an available field capacity larger than 20 mm, soil surface resistance is generally 100 s/m. For soils with smaller available field capacities, soil surface resistance decreases with increasing soil moisture, reaching 100 s/m at complete saturation (slightly smaller values are possible for soil moisture slightly below saturation).

import numpy
from matplotlib import pyplot
from hydpy.models.evap import *

pyplot.figure(figsize=[6.4, 3.2]).set_layout_engine("tight")
parameterstep()
nmbhru(1000)
soil(True)
wiltingpoint(20.0)
factors.soilwater = numpy.linspace(15.0, 45.0, nmbhru.value)
for fc in [20.0, 30.0, 40.0, 41.0]:
    fieldcapacity(fc)
    model.calc_soilsurfaceresistance_v1()
    pyplot.plot(
        factors.soilwater.values,
        factors.soilsurfaceresistance.values,
        label=f"FC-WP = {int(fc - wiltingpoint.values[0])} mm",
    )
pyplot.legend(frameon=False)
pyplot.xlabel("SoilWater [mm]")
pyplot.ylabel("SoilSurfaceResistance [s/m]")

That means "normal" bare soils always evaporate as much water as shallow bare soils that are entirely wet. Is this another possible cause for high evapotranspiration estimates?

[tyralla]

aerodynamic resistance

Besides the roughness length issue, the calculation of the aerodynamic resistance includes an unfortunate discontinuity regarding crop height (already mentioned in the documentation, but a figure might make that clearer):

import numpy
from matplotlib import pyplot
from hydpy.models.evap import *

pyplot.figure(figsize=[6.4, 3.2]).set_layout_engine("tight")
from hydpy.models.evap import *
parameterstep()
nmbhru(1)
hrutype(0)
derived.moy(0)
factors.windspeed10m = 3.0
cropheights = numpy.linspace(0.0, 15.0, 1000)
aerodynamicresistances = []
for ch in cropheights:
    cropheight.values = ch
    model.calc_aerodynamicresistance_v1()
    aerodynamicresistances.append(factors.aerodynamicresistance.values[0])
pyplot.plot(cropheights, aerodynamicresistances)
pyplot.xlabel("CropHeight [m]")
pyplot.ylabel("AerodynamicResistance [s/m]")

[tyralla]

PWP vs. PY

There is some confusion concerning the current lland parameters PWP and PY. We currently use PWP for calculating the soil surface resistance and PY for calculating the surface resistance of vegetation, water, and sealed surfaces. Besides that, lland uses PWP for calculating interflow and base flow.

First, we already agreed to use the same parameter to calculate both resistance terms. Hence, our new AET submodel requires only one parameter. I suggest naming it WiltingPoint. Unfortunately, this changes the old model behaviour for different PWP and PY settings.

Second, we must decide if PWP (responsible for runoff generation) and PY (responsible for evapotranspiration) address the same physical property. If so, lland_v3 and lland_v4 could automatically set the AET submodel parameter WiltingPoint based on previously defined PWP values. (Model users could still override this behaviour by defining different WiltingPoint values within control files.)

[tyralla]

"landuse" surface resistance

Maybe, the idea of understanding PY as the permanent wilting point is not so good. The following graph shows that the surface resistance of vegetation initially only increases slowly when the soil moisture starts to sink below this threshold (100 mm).

import numpy
from matplotlib import pyplot
from hydpy.models.evap import *

pyplot.figure(figsize=[6.4, 3.2]).set_layout_engine("tight")
parameterstep()
nmbhru(1000)
hrutype(0)
soil(True)
conifer(False)
derived.moy(0)
surfaceresistance(100.0)
wiltingpoint(100.0)
fieldcapacity(200.0)
factors.soilwater = numpy.linspace(0.0, 150.0, nmbhru.value)
model.calc_landusesurfaceresistance_v1()
pyplot.plot(factors.soilwater.values, factors.landusesurfaceresistance.values)
pyplot.ylim((0.0, 1000.0))
pyplot.xlabel("SoilWater [-]")
pyplot.ylabel("LanduseSurfaceResistance [s/m]")

There is another complication for coniferous forests (not for mixed forests). Only for this kind of vegetation "suboptimal" air temperature and humidity conditions can further reduce surface resistance. Generally, low air humidity results in increased surface resistance. For humid air, decreasing air temperature also results in increased surface resistance. For an infinite air temperature and air humidity at saturation, coniferous forests behave like, for example, deciduous forests.

import numpy
from matplotlib import pyplot
from hydpy.models.evap import *

pyplot.figure(figsize=[6.4, 3.2]).set_layout_engine("tight")
parameterstep()
nmbhru(100)
hrutype(0)
soil(True)
conifer(True)
derived.moy(0)
surfaceresistance(100.0)
wiltingpoint(100.0)
factors.soilwater = 100.0
humidities = numpy.linspace(0.0, 1.0, nmbhru.value)
for t in [-5.0, 0.0, 5.0, 10.0, 15.0, 20.0, 25.0, 30.0]:
    factors.airtemperature.values = t
    model.calc_saturationvapourpressure_v1()
    factors.actualvapourpressure = humidities * factors.saturationvapourpressure
    model.calc_landusesurfaceresistance_v1()
    pyplot.plot(
        humidities, factors.landusesurfaceresistance.values, label=f"T = {int(t)} °C"
    )
pyplot.legend(frameon=False)
pyplot.xlabel("relative humidity [-]")
pyplot.ylabel("LanduseSurfaceResistance [s/m]")

[tyralla]

leaf area index

lland_v3 and lland_v4 require the leaf area index for their AET and snow methods. After extracting both processes into submodels, there is no chance of sharing this information with the current HydPy mechanisms (except for the second-best solution of defining LAI values in auxiliary files).

I suggest we ignore this problem for now and duplicate the leaf area index parameter for now. Information sharing would not be straightforward even between a main model and its submodel for this case because it depends on the month and the land type (other models might favour other dependencies). So we better solve this problem later in a separate attempt.

[tyralla]

actual surface resistance

lland_v3 and lland_v4 determine a weighted average of soil and land use surface resistance that depends on the leaf area index and the possible sunshine duration. HydPy's and LARSIM's documentations reference Grant (1975), but I see only a loose connection at first glance. The following graph illustrates these dependencies of a simulation step of 24 hours. Still, it would be identical for any other time step (12 hours possible sunshine duration for a daily simulation is handled precisely like 30 minutes possible sunshine duration for an hourly simulation).

import numpy
from matplotlib import pyplot
from hydpy.models.evap import *

pyplot.figure(figsize=[6.4, 3.2]).set_layout_engine("tight")
parameterstep()
nmbhru(1)
hrutype(0)
soil(True)
derived.moy(0)
derived.hours(24)
factors.landusesurfaceresistance = 100.0
factors.soilsurfaceresistance = 200.0
leafareaindices = numpy.linspace(0.0, 10.0, 100)
for psd in [0.0, 6.0, 12.0, 18.0, 24.0]:
    inputs.possiblesunshineduration = psd
    actualsurfaceresistances = []
    for lai in leafareaindices:
        leafareaindex.values = lai
        model.calc_actualsurfaceresistance_v1()
        actualsurfaceresistances.append(factors.actualsurfaceresistance.values[0])
    pyplot.plot(leafareaindices, actualsurfaceresistances, label=f"PSD = {int(psd)} h")
pyplot.legend(frameon=False)
pyplot.xlabel("LeafAreaIndex [-]")
pyplot.ylabel("ActualSurfaceResistance [s/m]")

[tyralla]

soil heat flux

There is currently some confusion in lland_v3 and lland_v4 because of two different definitions of soil heat fluxes. G is related to the Penman-Monteith method, and WG is related to the Knauf method. Extracting the Penman-Monteith-related equations simplifies lland because we can remove G.

However, evap also defines a corresponding parameter called SoilHeatFlux. Unfortunately, positive G values increase evaporation estimates (in lland) while positive SoilHeatFlux values decrease evaporation estimates (in evap). The SoilHeatFlux definition seems more common (and agrees with LARSIM's documentation). If I member correctly, we decided on the definition of G for consistency with the definition of WG. As this argument becomes irrelevant, I decided to stick to the SoilHeatFlux definition and adjust the methods to move from lland to evap accordingly.

[tyralla]

Penman vs. Penman-Monteith

lland_v3 and lland_v4 use Penman (always based on daily input data) for estimating evaporation from water areas, Penman-Monteith with zero surface resistance for interception evaporation and Penman-Monteith with properly calculated surface resistance for evapotranspiration from soils. In the documentation on the methods Return_Penman_V1 and Return_PenmanMonteith_V1, we compare both approaches for similar (open water surface) conditions. For pure energy forcings, Return_Penman_V1 produces slightly higher estimates. For pure wind forcings, Return_PenmanMonteith_V1 produces much higher estimates. The only error in the comparison could be inconsistency between the wind speed (used by Return_Penman_V1) and the aerodynamic resistance (used by Return_PenmanMonteith_V1). But, maybe, this could be another cause for overestimation tendencies.

Regarding the aerodynamic resistance: I did not check, but I guess I took the wind speed at 2 m used by Penman, adjusted it to 10 m with Calc_WindSpeed10m_V1, and finally applied Calc_AerodynamicResistance_V1 assuming a crop height of zero.

[tyralla]

Wigmosta et al. (1994)

lland_v3 and lland_v4 use the following equation, based on Wigmosta et al. (1994), to reduce soil evapotranspiration in case of simultaneous interception evaporation:

$SoilEvapotranspiration{new} = \frac{PotentialInterceptionEvaporation - InterceptionEvaporation} {PotentialInterceptionEvaporation} \cdot SoilEvapotranspiration{old}$

It is unclear how to deal with (partly) negative values. So far, the equation is applied without modification, except for zero potential interception evaporation values, for which soil evaporation is also set to zero.

[tyralla]

soil heat flux

lland_v3 and lland_v4 calculate two soil heat fluxes. One is MORECS-like and targets evapotranspiration, while the other is Knauf-like and targets the energy balance of snow layers. The separation of these two fluxes into two models is an immense simplification. However, for now, I will still use the workaround method Calc_G_V2 instead of the problematic method Calc_G_V1.

We should come back to this later. Here are only some additional remarks:

• The MORECS report states that the night-time soil heat flux calculation relies on a scheme suggested by Wales-Smith and Arnott (1980). This paper would possibly help, but it seems to be unpublished.
• The mentioned approach relies on monthly values of "average daily heat storage in soil" the authors calculated based on soil temperature measurements and estimates of soil thermal capacities at Kew (London).
• For vegetation other than grass, the day-time soil heat flux depends on the leaf area index: $Gd = (0.3 - 0.03 \cdot LAI) \cdot R{nd}$ So the daily soil heat flux could be zero or even negative for dense forests?

[tyralla]

daily averages

Here, I assume a simulation step size of one hour for simplicity.

lland_v3 and lland_v4 apply the Penman equation based on the aggregated (averaged or summed) input data of the last 24 hours. Additionally, they insert the "daily" net longwave radiation into the Penman-Monteith equation, which they calculate based on other aggregated input data (air temperature, actual vapour pressure, sunshine duration, and possible sunshine duration).

Some notes:

• Calc_TemLTag_V1 averages the subbasin-wide, uncorrected air temperature and Calc_TKorTag_V1 adjusts this daily average to the individual zone. In evap_morecs, we must average the already adjusted zone-specific air temperatures instead because evap_morecs does not know TemL (and we do not want to lose chances of coupling evap_morecs with other models that perform different temperature adjustments). However, due to the simple, additional temperature adjustment, we still keep compatibility with the old lland_v3 and lland_v4 results.
• lland_v3, lland_v4, and evap_morsim aggregate the following properties in the same by calculating their mean or sum values over the last 24 hours: sunshine duration, possible sunshine duration, wind speed at 2 m above ground, and soon relative humidity (see the following note).
• For the actual vapour pressure, the saturation vapour pressure, and the saturation vapour pressure slope, lland_v3 and lland_v4 follow another strategy. They first calculate the daily saturation vapour pressure and the daily saturation vapour pressure slope directly based on the daily adjusted air temperature. Then, they average the relative humidity values of the last 24 hours. Finally, they use the derived daily relative humidity and saturation vapour pressure values to calculate the daily actual vapour pressure. The same strategy should work for evap_morecs as well. I will try it.
• As explained above, the daily net longwave radiation is always calculated directly. Neither the Penman nor the Penman-Monteith equations rely on hourly net longwave radiation values (this is different for the snow module, but that is another topic).
  • Daily net shortwave radiation depends on daily global radiation, which is also determined by averaging hourly values. Confusingly, it also depends on ActualAlbedo and not something like DailyActualAlbedo (which does not exist). On first thought, I suspected an unnecessary inconsistency here. On second thought, I am unsure because the switching between snow and "normal" evaporation calculations complicates matters. The general handling of snow-related albedo modification is still undecided, so I leave this problem for later.

[tyralla]

actual albedo - part 2

I decided to introduce a SnowAlbedoModel_V1 submodel interface. lland_v3 and lland_v4 comply with this interface. Also, lland_v3 and lland_v4 now calculate the actual albedo only if zones are covered with snow. Otherwise, they set it to nan.

The new thing regarding our current submodel strategy is that the new get_snowalbedo method of lland_v3 and lland_v4 sometimes returns nan and other times relevant values. Hence, evap_morsim checks the given value and uses it if it is not nan. But if it is nan, indicating snow-free conditions, it takes the land type- and month-specific value from the Albedo parameter.

In this regard, evap_morsim considers the snow albedo submodel already as optional. I stretched this further by officially defining it as optional. So, models like hland_v1 or lland_v1, which do not model the snow albedo, should also be able to use evap_morsim. For such combinations, evap_morsim always uses albedo for snow-free surfaces. (Except one would add a real submodel that provides independent snow albedo estimates.)

(Theoretically, evap_morsim could apply a fixed snow albedo value if hland_v1 or lland_v1 say that there is a snow cover via the SnowCoverModel_V1 interface.)

For now, this mechanism should work. However, we probably will move the rest of the actual albedo calculations to a new snow submodel as soon as we start extracting snow processes from our main models. Then, the "use the main model as a submodel" mechanism will not apply anymore.

[tyralla]

interception condensation

The Penman-Monteith equation can return negative interception evaporation values. The corresponding increase in the amount of intercepted water can result in an exceedance of the interception capacity. Currently, lland_v3 and lland_v4 accept such a threshold violation for the simulation time step where it initially occurs. They only convert this excess to stand precipitation at the beginning of the next time step.

This behaviour differs from the new behaviour of lland_v1 and the hland application models, which directly convert the excess condensation to stand precipitation, which seems more logical. So I intended to adjust lland_v3 and lland_v4 accordingly. However, that evap_morecs does not (always) simulate interception evaporation if there is a snow cover complicates the issue because a snow cover could originate from condensation-based throughfall.

For clarification, here are two possible workflows.

1. Allow excess:

   • Convert precipitation to throughfall
   • Add the (frozen part of the) throughfall to the snow layer
   • Calculate interception evaporation/condensation (if not suppressed by snow generated at this time step or earlier)
   • Do not modify too-high amounts of intercepted water

2. Prohibit excess:

   • Convert precipitation to throughfall
   • Calculate interception evaporation/condensation (if not suppressed by snow generated earlier)
   • Add any condensation excess to the previously generated throughfall
   • Add the (frozen part of the) throughfall to the snow layer

I favour the second approach, as it cannot cause writing inconsistent initial condition files. From other viewpoints, both have inconsistencies, which likely cannot be avoided when using such an ad-hoc integration strategy. However, at least for hourly time steps, the effects of these inconsistencies should be of minor importance.

Maybe it seems contradictory, but I will still execute the soil evapotranspiration routine after performing the snow routine for now so that a freshly built-up snow layer can directly suppress evapotranspiration from the soil. I am afraid we need to revisit our decisions regarding execution orders anyhow, so this topic will likely come up automatically.

[tyralla]

interception condensation / snow albedo / snow cover

Calculating interception evaporation and condensation before executing the snow module raises another problem: evap_morsim then queries the current snow albedo, which is still not determined by lland_v3 or lland_v3 within the current time step. Hence, evap_morsim usually receives the value calculated in the last step. However, such behaviour is unacceptable as the value of the first simulation step would depend on previous work done with the model.

This problem is currently not relevant for lland_v1 and the hland models as they do not provide snow albedo values. However, it might become relevant when extracting the snow routines so that they could use Knauf's snow energy balance method or something similar.

While thinking, I realise the same problem exists for the current degree of snow cover, which evap_morsim queries from lland_v3 or lland_v4 to suppress interception evaporation for non-forest land types eventually. So, unfortunately, coupling evap_morsim to lland_v1 or a hland model comes with the same problem.

I am afraid there is no chance of a completely satisfactory solution. I see two options:

• Let all hland and lland models first calculate the snow routine and then the evaporation routines, including interception evaporation.
• Define two distinct interfaces. One expects the order: interception evaporation -> snow routines -> soil evapotranspiration. The other expects the order: snow routines -> interception evaporation -> soil evapotranspiration.

On first thought, I prefer the first option because it appears less confusing. But we would need to adjust lland_v1 and the hland models regarding interception condensation. Then, we could either disallow any condensation (which does not seem right to me) or always leave the condensed water for the next time step (which means we could not trim the interception storages regarding their maximum values anymore, but regarding simulations, this change should only show minimal effects).

[tyralla]

interception condensation / snow albedo / snow cover - part 2

While thinking, I realise the same problem exists for the current degree of snow cover, which evap_morsim queries from lland_v3 or lland_v4 to suppress interception evaporation for non-forest land types eventually. So, unfortunately, coupling evap_morsim to lland_v1 or a hland model comes with the same problem.

That is only half-true. As long as the hland and lland models calculate the snow cover degree based on the current snow states on the fly when get_snowcover is called (which is currently the case), the queried snow cover is always adequately defined.

However, the reasoning above remains the same. Someday, one might turn "snow cover" into a factor sequence, which then would actually introduce the problem - possibly for a long time until someone notices. So, trying out a fixed "first snow, then evaporation" order still seems preferable.

[tyralla]

to do

• [x] Assume that the field capacity of lland_v3 and lland_v4 is the same as that of evap_morsim?
• [x] Assume that the permanent wilting point of lland_v3 and lland_v4 is the same thing as the wilting point (or PY) of evap_morsim?
• [x] Let evap_morsim and its main model share their standard scalar input time series (sunshine duration, wind speed, relative humidity, possible sunshine duration, global radiation)? (Maybe we should postpone this decision because the snow modularisation will remove many input time series from lland, and it is currently unclear what turning the meteo models into submodels will bring.)
• [x] Now that we could couple, for example, lland_v3 with evap_minhas, we should probably add the same "temporal smoothing mechanism" as for lland_v1 for preventing the possibility of strong diurnal patterns in simulated discharge. Or should we move the mechanism to evap_minhas and evap_aet_hbv96?
• [x] evap_morsim uses log sequences for calculating "daily" averages or sums (for the last 24 hours). Hence, using evap_morsim operationally requires that we consider submodels' log (and state) sequences when reading and writing condition files, similar to reading and writing control files.
• [x] Decide whether snow interception should suppress evapotranspiration from the soil (see the snow-in-tree-canopies example). If not, we should explain why or at least remove the "snow-in-tree-canopies" sentence.
• [x] Let lland and evap automatically share information on the measuring height of wind speed?
• [x] Let lland and evap automatically share information on the month- and land-type specific leaf area index?

[tyralla]

• Let evap_morsim and its main model share their standard scalar input time series (sunshine duration, wind speed, relative humidity, possible sunshine duration, global radiation)? (Maybe we should postpone this decision because the snow modularisation will remove many input time series from lland, and it is currently unclear what turning the meteo models into submodels will bring.)

This should be sufficiently solved. See issues #75 and #124.

[tyralla]

• Assume that the field capacity of lland_v3 and lland_v4 is the same as that of evap_morsim?
• Assume that the permanent wilting point of lland_v3 and lland_v4 is the same thing as the wilting point (or PY) of evap_morsim?

We agreed that PY is different from the wilting point. We now use the SoilMoistureLimit parameter to represent MORECS' PY parameter in evap_morsim. Note that we also use SoilMoistureLimit to represent HBV's LP parameter in evap_aet_hbv96.

Correspondingly, we use the more neutral parameter MaxSoilWater instead of FieldCapacity, which is the same as HydPy-L's WMax parameter. Hence, lland_v3 and lland_4 should be able to pass this information to evap_morsim automatically.

The equations become

$$ SoilSurfaceResistance = \begin{cases} 100 &|\ W_m> 20 \ \frac{100 \cdot W_m}{W_a + W_m/100} &|\ W_m\leq 20 \end{cases} $$

$$ W_a = SoilWater $$

$$ W_m = MaxSoilWater $$

and

$$ LanduseSurfaceResistance = r^* \cdot \begin{cases} 3.5 \cdot (1 - W / \tau) + exp(0.2 \cdot \tau / W)&|\ W \leq \tau \ exp(0.2) &|\ W > \tau \end{cases} $$

$$ r^*= \begin{cases} r&|\ \overline{C} \ 10\,000 &|\ C \ \land \ (T \leq -5 \ \lor \ \Delta \geq 20) \ min\left(\frac{25 \cdot r}{(T + 5) \cdot (1 - \Delta / 20)}, \ 10\,000\right) &|\ C \ \land \ (-5 < T < 20 \ \land \ \Delta < 20) \ min\left(\frac{r}{1 - \Delta / 20}, \ 10\,000\right) &|\ C \ \land \ (20 \leq T \ \land \ \Delta < 20) \end{cases} $$

$$ W = SoilWater $$

$$ C = Conifer $$

$$ \tau = SoilMoistureLimit \cdot MaxSoilWater $$

$$ r = SurfaceResistance $$

$$ \Delta = SaturationVapourPressure - ActualVapourPressure $$

$$ T = AirTemperature $$

[tyralla]

• Decide whether snow interception should suppress evapotranspiration from the soil (see the snow-in-tree-canopies example). If not, we should explain why or at least remove the "snow-in-tree-canopies" sentence.

The current behaviour seems inconsistent but is the best agreement with LARSIM we could achieve. If no one objects, I will leave it like this and just remove the to-do sentence. (discussed offline with @sivogel)

[tyralla]

• Let lland and evap automatically share information on the measuring height of wind speed?
• Let lland and evap automatically share information on the month- and land-type specific leaf area index?

This would remove a lot of redundancies in the affected control files, especially regarding the leaf area index. So we decided on a Yes.

[tyralla]

• Now that we could couple, for example, lland_v3 with evap_minhas, we should probably add the same "temporal smoothing mechanism" as for lland_v1 for preventing the possibility of strong diurnal patterns in simulated discharge. Or should we move the mechanism to evap_minhas and evap_aet_hbv96?

We thought about this thoroughly and moved the evaporation-damping mechanism into evap_m and evap_mlc, which convert (grass) reference evapotranspiration to potential evapotranspiration. This is the suitable place for adjusting "land evapotranspiration" to "water evaporation" and gives us the most flexibility to apply this damping when using different model combinations.

[tyralla]

So, everything seems settled. I will run the tests a last time and push the changes into the master branch after that.