search

hydpy-dev / hydpy

A framework for the development and application of hydrological models based on Python
GNU Lesser General Public License v3.0
35 stars 17 forks source link

Improve HydPy-Evap-MORSIM or implement HydPy-Evap-AMBAV-1.0 #118

Closed tyralla closed 6 months ago

tyralla commented 12 months ago

We implemented HydPy-Evap-MORSIM as discussed in #111. There is still room for minor technical improvements, but the model's "hydrological functioning" is clear now. We see a marked overprediction in evaporation in a project area in northern Germany, which agrees with LARSIM's tendency to estimate high (potential) evapotranspiration rates. Considering also the relevant uncertainties in some process equations (partly discussed in #111), we decided to be inspired by the original AMBAV model, developed and used by the German Meteorological Service (Deutscher Wetterdienst). For a start, we will create a new model called HydPy-Evap-AMBAV-1.0 or evap_ambav_1, which works precisely like HydPy-Evap-MORSIM, and then replace some of its equations with more AMBAV-like equations. I will note the changes in this issue to keep the overview and everybody informed. We will see if this effort results in a complete AMBAV-like model or more in a MORSIM-AMBAV hybrid.

tyralla commented 12 months ago

This is the list of MORSIM equations we need to remove, modify, or replace with other approaches eventually (I add no comments if we can leave it like it is):

tyralla commented 12 months ago

Soil heat flux

This is a good starting point as method Calc_SoilHeatFlux_V3 is just a stopgap.

The AMBAV equation is based on the previously calculated net radiation ($R_n$) and is simple:

$$ G = (a - 0.03 \cdot LAI) \cdot R_n $$

However, AMBAV's documentation suggests at least defining different values of $a$ during daytime (possibly 0.2) and nighttime (0.5). So, we should allow setting two different values for day and night.

The AMBAV documentation also mentions that $a$ is 0.0 during the first two hours after sunrise and often 0.8 for wet soils. The first thought is surprising, and the second is too vague to be put into an equation without further assumptions. Hence, I will ignore these two statements for now.

AMBAV seems to calculate the time of sunrise ($SA$) and sunset ($SU$) for distinguishing between day and night with the following equations:

$$ SA = \frac{12}{\pi} \cdot \arccos \left( tan \ \theta \cdot \tan \ \phi + \frac{0.0145}{cos \ \theta \cdot cos \ \phi}\right) $$

$$ SU = 24 - SA $$

This is nearly the same as the FAO56 and already implemented by method Calc_TimeOfSunrise_TimeOfSunset_V1. Unfortunately, this equation is part of HydPy-Meteo and not of HydPy-Evap, implemented by the FAO56-compliant models meteo_v003 and meteo_v004. I prefer to avoid copy-pasting these astronomical calculations into HydPy-Evap, meaning HydPy-Evap-AMBAV-1.0 relies on one of these submodels or external sunrise and sunset estimates. (Model configuration will become complicated, but things will significantly improve as soon as we convert the stand-alone meteo models to submodels.)

tyralla commented 12 months ago

Soil heat flux, part 2

The equation

$$ G = (a - 0.03 \cdot LAI) \cdot R_n $$

should also have been included in Evap-MORSIM via method Calc_G_V1. However, Calc_G_V1 suffers two or three problems, which is why we used Calc_G_V2 instead.

One of these occurs for vegetation with large leaf area indices. MORSIM would have used $a=0.3$, and AMBAV seems to use $a=0.2$. The leaf area index can be larger than ten for tree-like vegetation so that the soil heat flux would be negative during the day. This may be fine for the original AMBAV model, as it focuses on crop plants like wheat and corn.

Theoretically, we could just prevent negative subterms:

$$ G = max(a - 0.03 \cdot LAI, \ 0) \cdot R_n $$

However, that would be inconsistent with the calculation for the nighttime, which is either

$$ G = 0.5 \cdot R_n $$

or

$$ G = (0.5 -0.03 \cdot LAI) \cdot R_n $$

(I think the latter).

$a=0.2$ is said to hold for grass. Maybe someone determined values suitable for woods?

If no one suggests a better solution, I would go for the following equations:

$$ G_{day} = max(0.2 - 0.03 \cdot LAI, \ 0) \cdot R_n $$

$$ G_{night} = max(0.5 - 0.03 \cdot LAI, \ 0) \cdot R_n $$

tyralla commented 12 months ago

Soil heat flux, part 3

Two updates.

First, the implementation is simpler than expected because we can directly use the astronomically possible sunshine duration instead of the times of sunrise and sunset, which is already known by HydPy-Evap-MORSIM.

Second, when writing the unit tests, I could not find a good excuse for why woods should suppress soil head fluxes entirely during the day but not at night. So, I modified the equation a little in a manner that seems more consistent to me:

$$ G = max(a - b \cdot LAI, \ 0) \cdot R_n $$

$$ a = \omega \cdot 0.2 + (1 - \omega) \cdot 0.5 $$

$$ b = a \cdot 0.03 / 0.2 $$

$$ G = SoilHeatFlux $$

$$ LAI = LeafAreaIndex $$

$$ R_n = NetRadiation $$

$$ \omega = PossibleSunshineDuration / Hours $$

tyralla commented 12 months ago

Wind speed / Aerodynamic resistance

Nearly twenty years have passed since I visited some introductory seminars in meteorology. I do my best to collect the relevant information, but help from a meteorologist would be highly appreciated...

First, AMBAV assumes neutral atmospheric conditions, so we can principally use the logarithmic wind profile implemented by method Return_AdjustedWindSpeed_V1.

One thing that confuses me is that the AMBAV documentation states we cannot neglect zero plane displacement ($d$) for woods and suggests using a simple correction following Oliver (1974):

$$ u(10 + d) = 0.6 \cdot u(10) $$

This is the only topic where the documentation mentions woods, so why here? And while asking myself if we should use it (I tend to assume so), I come to the second topic.

This is how AMBAV calculates the aerodynamic resistance:

$$ r_a = \frac{l}{\kappa^2 \cdot u(z)} \cdot ln \ \left( \frac{z - d}{z_0} \right) + ln \ \left( \frac{z - d}{(z0){H,E}} \right) $$

$r_a$: aerodynamic resistance (really, the documentation uses $(ra){H,E}$ instead of $r_a$ here, which is used in the Penman-Monteith equation. Do both terms mean the same here?)

$l$: ???

$\kappa$: Karman constant (should we use 0.4?)

$z$: height (for us, always two meters?)

$d$: zero plane displacement (see below)

$z_0$: roughness length for momentum exchange (see below)

$(z0){H,E}$: roughness length for sensible and latent heat exchange ($0.2 \cdot z_0$)

Furthermore, the AMBAV documentation mentions multiple equations for calculating the roughness length and zero plane displacement based on the crop height and/or leaf area index. Which should we choose?

The documentation nowhere mentions $l$. Possibly, I could solve the mystery by digging into the equations, but I am thankful for any hint...

tyralla commented 12 months ago

Albedo

AMBAV calculates the current total albedo ($\alpha$) by weighting the soil's albedo ($\alpha_s$) and the plant's albedo ($\alpha_p$):

$$ \alpha = \alpha_s + 0.25 \cdot (\alpha_p - \alpha_s) \cdot LAI \ | \ LAI < 4
$$

$$ \alpha = \alpha_p \ | \ LAI \ge 4 $$

For crops, $\alpha_p$ is said to only vary between 0.21 and 0.26 (in Germany), but $\alpha_s$ receives special attention.

First, it is said to vary between 0.15 for dark soils and 0.38 for light sands.

Second, it is corrected to 0.8 for complete snow covers and $(\alpha_s + 0.8) / 2$ for incomplete snow covers. (Why does snow not affect the plant albedo?)

Third, it is halved for wet soils. A soil is said to be wet if the rainfall of the last 24 hours was at least 0.5 mm, 1.0 mm, or 1.5 mm for different times of the year.

The documented strategies differ from how things are implemented in Evap-MORSIM, which, for example, deals with fixed monthly total albedo values. Here, we would need to introduce two new parameters, AlbedoSoil and AlbedoPlant. The values of AlbedoPlant could be set based on the defined land use types (HRUType). AlbedoSoil could be selected based on a new SoilType parameter (Should we then rename HRUType to LandType for consistency? Probably depends on whether we accept WATER as a type of land...). Additionally, we would require a seasonal threshold parameter for distinguishing between wet and dry soils.

tyralla commented 12 months ago

Longwave radiation

Regarding surface radiation, the AMBAV documentation suffers from some confusing references. It says AMBAV relies on equation 2.5.17, which is not defined. I assume equation 2.5.18 is meant ($AirTemperature$ in Kelvin):

$$ RL_{s} = 0.97 \cdot \sigma \cdot (AirTemperature + 293.1)^4 + 0.07 \cdot (1 - \alpha) \cdot GlobalRadiation $$

The equation for the atmospheric counter radiation cannot be correct:

$$ RL_{a} = (1 + C \cdot N^2) \cdot \big(1 - 0.261 \cdot exp(-7.77 \cdot 10^{-4} + AirTemperature^2)\big) \cdot \alpha \cdot (AirTemperature + 273.1)^4 $$

I found a comparison of different equations for atmospheric emissivity (table 1). It seems we must simply turn a $+$ into a $\cdot$:

$$ RL_{a} = (1 + C \cdot N^2) \cdot \big(1 - 0.261 \cdot exp(-7.77 \cdot 10^{-4} \cdot AirTemperature^2)\big) \cdot \alpha \cdot (AirTemperature + 273.1)^4 $$

We can implement both equations as noted. However, $N$, the sky coverage fraction, and $C$, a coefficient reflecting the cloud type, are often not directly available for simulation.

Regarding $N$, we could either use the relation between global radiation and clear sky solar radiation (as in Calc_NetLongwaveRadiation_V1) or sunshine duration and astronomically possible sunshine duration (as in Calc_RLAtm_V1). Both work only during daytime, of course. Which approach should we favour, which equation should we use, and how should we estimate the nighttime values?

Regarding $C$, I am not aware of an established estimation procedure.

tyralla commented 11 months ago

Surface resistance

A little surprise: AMBAV does not consider water deficits when determining the surface resistance. In contrast to Evap-MORSIM, it first calculates potential evapotranspiration and reduces it to actual evapotranspiration in a subsequent step (see below). Hence, it appears reasonable to implement two separate application models, HydPy-Evap-PET-AMBAV-1 and HydPy-Evap-AET-AMBAV-1 as we did for HydPy-Evap-PET-HBV96 and HydPy-Evap-AET-HBV96.

Like Evap-MORSIM, AMBAV calculates the total surface resistance ($rs$) based on Grant (1975) by weighting the plant's surface resistance ($rs_p$) and the soil's surface resistance ($rs_s$):

$$ \frac{l}{rs} = \frac{l-I}{rs_p} + \frac{I}{rs_s} $$

However, $I$ is not always calculated as $0.7^{LAI}$ as in Evap-MORSIM but as $0.8^{LAI}$ for leaf area index values smaller than one.

Here, the mysterious $l$ emerges again. After going through the complete documentation, I think it is just a formatting problem and $l=1$, so that:

$$ \frac{1}{rs} = \frac{1-I}{rs_p} + \frac{I}{rs_s} $$

During the nighttime, $rs_p$ is 2800 s/m. We could weight the day and nighttime values of simulation steps around sunrise and sunset as we do in Calc_LanduseSurfaceResistance_V1.

The calculation of $rs_s$ relies on the same season-specific precipitation thresholds as the calculation of $\alpha_s$. No equation is given; I interpret the text as follows:

$$ rss^t = \begin{cases} rs{s,wet} &|\ Precipitation_t \geq Threshold \ rs_s^{t-1} + c &|\ Precipitation_t < Threshold \end{cases} $$

If this is what we want to do, what are suitable names and values for $rs_{s,wet}$ (WetSoilSurfaceResistance, 100 s/m?) and $c$ (SoilSurfaceResistanceIncrease, 10 s/m/day?)?

tyralla commented 11 months ago

Interception

The description of interception evaporation is relatively brief. Obviously, non-wetted surface fractions (where AMBAV applies $rs$ as explained above) and wetted surface fractions ($rs=0$) can coexist. The calculation of these fractions is less clear. The following is mostly guesswork:

$$ MaxInterceptedWater = 0.2 \cdot LeafAreaIndex $$

$$ WettedSurfaceFraction = w = InterceptedWater / MaxInterceptedWater = i / m $$

$$ PotentialEvapotranspiration = PET = w \cdot PET(0) + (1 - w) \cdot PET(rs) = w \cdot p_0 + (1 - w) \cdot p_r $$

So, instead of assuming the intercepted water is always distributed as uniformly as possible, we (and eventually AMBAV) would say that it is distributed as ununiformly as possible.

If we decide on this approach, the potential evaporation rate would depend on the current amount of intercepted water. We could probably solve the corresponding ODE without errors based on the following (still unchecked) equation:

$$ i_t = c \cdot exp(t \cdot (p_r - p_0) / m) - (p_r \cdot m)/(p_s - p_r) $$

The above solution assumes that $w \cdot p_0 + (1 - w) \cdot p_r$ is subtracted from $i$, which is confusing. However, subtracting only $w \cdot p_0$ from $i$ would imply $i$ reaches zero only asymptotically if we solve the ODE without errors.

tyralla commented 11 months ago

Reduction of potential to actual evapotranspiration

AMBAV's approach of reducing potential evapotranspiration to actual evapotranspiration is very unique. Basically, it follows Slabbers (1980). He calculates $f$, the fraction of available soil moisture at which the reduction starts, via:

$$ f = 0.94 +0.26 \cdot \frac{\psi_p^{crit}}{PET} $$

$\psi_p^{crit}$ is the critical leaf water potential (beware of the inconsistencies in its signum!).

Principally, AMBAV applies $f$ by comparing it with the fraction of the currently available water content ($AWC$) and the highest possible available water content ($MaxAWC$):

$$ AET = \frac{AWC}{min(f \cdot MaxAWC, \ AWC)} \cdot PET $$

However, the AMBAV documentation states that this approach would result in zero reductions if PET is below 4 mm/d for some stands, so the soil's available water content could become negative. AMBAV's solution seems to be to recalculate $f$ based on another equation if the original equation suggests that a reduction should only take place for $AWC < 0.4 \cdot MaxAWC$:

$$ f{mod} = \begin{cases} f{orig} &|\ f_{orig} \geq 0.4\ 0.727 \cdot \frac{PET}{\psip^{crit} + 0.05} &|\ f{orig} < 0.4 \end{cases} $$

Slabbers's method targets daily simulation steps. Can we adjust it to sub-daily (e.g. hourly) intervals?

We would need to introduce $\psi_p^{crit}$ as a now land type-specific (and eventually also month-specific?) parameter. Are suitable aggregated values for land cover types such as "field" or "coniferous forest" available?

Does this approach work for tree-like vegetation?

tyralla commented 11 months ago

Root depth

AMBAV 2.0 separates the soil into multiple layers, but AMBAV 1.0 seems to follow a more lumped approach. At least conceptionally, it separates the soil into two compartments. The first corresponds to the rooted depth, and the second lies below (but at a maximum depth of 60 cm?). Water is vertically evenly distributed within the rooted depth, and differences in root density are ignored. Precipitation only adds to this first compartment until it reaches the available water capacity. After that, it immediately percolates to the second compartment.

Eventually, the second compartment is separated into layers. Then, the first layer would only route water if its own field capacity is exceeded. However, I think AMBAV 1.0 does not implement such layers.

Following this concept, plants have little water supply when they start to grow and successively make (already available) soil water accessible for evapotranspiration during their growth. This is said to prevent overestimating evapotranspiration in April, May, and, for some plants, June (in Germany).

This is probably the only idea of AMBAV, whose implementation into Evap-AMBAV would be a little more complicated. We would need to think about how water fluxes like interflow or base flow (calculated by main models like HydPy-L-Land) should modify the different compartment contents.

Technically, this is similar to and simpler than the coupling between HydPy-L and GARTO we already implemented (#89). However, we first need to talk about whether we require such a separation into two compartments at all and whether we should make their size depend on the time-variable rooted depth.

So far, no model implemented in HydPy considers root growth. In case we want to change this: Is it okay to implement it with a month- and land use-specific parameter (like for the leaf area index), or could this result in sudden jumps in evapotranspiration?

tyralla commented 11 months ago

Water areas

Besides tree-like vegetation, AMBAV does not deal with water areas. Do we want to apply the Penman-Monteith equation with zero surface resistance values in such cases? Or should we use the original Penman equation?

Either way, we need to consider thermal inertia here.

Eventually, we can reuse most Penman-related methods already implemented for Evap-MORSIM.

tyralla commented 11 months ago

Wind speed / Aerodynamic resistance, part 2

Two conclusions and one additional question from a discussion we just had:

tyralla commented 11 months ago

Wind speed / Aerodynamic resistance, part 2

Our decisions:

tyralla commented 11 months ago

Wind speed / Aerodynamic resistance, part 3

I suggest to introduce a minimum value for the roughness length as follows:

$$ z_0 = max(0.13 \cdot h, \ 0.00013) $$

Such a boundary is required to avoid zero divisions when applying the logarithmic wind profile equation. 0.00013 lies between the values of water areas and bare soils given in the AMBAV documentation.

Should we introduce a maximum value as well?

tyralla commented 11 months ago

Wind speed / Aerodynamic resistance, part 4

To give an impression of our wind speed adjustment assumptions: With a wind speed of 5 m/s measured 10 m above a grass surface, we get a wind speed 2 m above ground (or zero displacement height) of 4.3 m/s for water, 4.0 m/s for grass (2.3 cm height), and 0.6 m/s for wood (10 m height).

tyralla commented 11 months ago

Wind speed / Aerodynamic resistance, part 5

In my understanding, we decided on the following equation for calculating the aerodynamic resistance:

$$ ra = \frac{1}{0.4^2 \cdot u{2m}} \cdot ln \ \left( \frac{z - d}{z_0} \right) + ln \ \left( \frac{z - d}{0.2 \cdot z_0} \right) $$

Can we assume that $z - d = 2m$? (at least, this would seem to be in line with how we handle forests):

$$ ra = \frac{1}{0.4^2 \cdot u{2m}} \cdot ln \ \left( \frac{2}{z_0} \right) + ln \ \left( \frac{10}{z_0} \right) $$

But, the AMBAV documentation defines $d$ values, so not using them seems wrong.

Additionally, for $z_0 \geq 2$, the first resistance term related to wind speed is zero. Hence, missing turbulence due to nearly zero wind speed would strongly reduce grass evapotranspiration but only slightly reduce forest evapotranspiration. Wind speed could even become entirely irrelevant for extremely high and rough canopies. Is this okay? Or should we at least introduce a maximum value for $z_0$, say 1.3 m?

Any suggestions are welcome...

tyralla commented 11 months ago

Wind speed / Aerodynamic resistance, part 6

The first results of calculating $r_a$ with the above equation were unconvincing. There is a typo in the AMBAV documentation's corresponding equation. Maybe I should interpret it so:

$$ ra = \frac{1}{0.4^2 \cdot u{2m}} \cdot \Bigg( ln \ \left( \frac{2}{z_0} \right) + ln \ \left( \frac{10}{z_0} \right) \Bigg) $$

At first sight, some results then look better. However, following Cleverly et al. (2013) it should be:

$$ ra = \frac{1}{0.4^2 \cdot u{2m}} \cdot \Bigg( ln \ \left( \frac{2}{z_0} \right) \cdot ln \ \left( \frac{10}{z_0} \right) \Bigg) $$

No fun here...

tyralla commented 11 months ago

Wind speed / Aerodynamic resistance, part 7

According to the FAO guidelines (neglecting zero plane displacement):

$$ r_a = \frac{ln(2 / z_0) \cdot ln(20 / z0)}{0.41^2 \cdot u{2m}} $$

(They say the first and second $z$ values are the measuring height of wind speed and humidity, but insert 2 meters in both cases.)

So, after rereading the AMBAV documentation, I have decided on the following equation:

$$ r_a = \frac{ln(2 / z_0) \cdot ln(10 / z0)}{0.4^2 \cdot u{2m}} $$

tyralla commented 11 months ago

Wind speed / Aerodynamic resistance, part 8

Adjusting the wind speed to ten meters is unnecessary and potentially harmful for large trees with rough canopies.

We better use:

$$ r_a = \left( \frac{ln(10 / z0)}{0.4} \right)^2 \cdot u{10m} $$

If I remember correctly, I found this equation also in Oke (1987).

I will use this equation. But then, I cannot reset low wind speeds at the height of 2 m to 0.5 m/s as easily as before. As it is still unclear whether we want to adopt this idea and for which height the minimum value of 0.5 m/s holds, I will skip this adjustment for now. (Such a minimum wind speed is also not mentioned in the AMBAV documentation.)

tyralla commented 11 months ago

Surface resistance, part 2

I introduce the parameters LeafResistance and DrySoilResistance. DrySoilResistance covers all hydrological response units from which water can evaporate directly from the soil, whereas LeafResistance only requires values for hydrological response units where plants grow. Both parameters do not allow to define different values for different seasons.

Currently, no flag is available to indicate "where plants grow", so I add the parameter Plant for this purpose (A flag that indicates whether the individual zones contain any vegetation.) and make sure HydPy-L(ARSIM) and HydPy-H(BV) set it automatically when using HydPy-Evap-AMBAV-1.

While introducing Plant, I found improving how main models and (sub-)submodels share information necessary, but I will discuss this in #90.

tyralla commented 11 months ago

Surface resistance, part 2

MORECS, LARSIM and Evap-MORSIM use the following equation for calculating the total surface resistance during daytime:

$$ \frac{1}{rs} = \frac{1-0.7^{LAI}}{rs_p} + \frac{0.7^{LAI}}{rs_s} $$

AMBAV relies on the same equation with a minor modification (discussed above).

However, I assumed AMBAV uses the same equation during nighttime, taking the very high plant-specific surface resistance of 2800 s/m due to closed stomata. However, MORECS, LARSIM, and Evap-MORSIM then use:

$$ \frac{1}{rs} = \frac{LAI}{2500} + \frac{1}{rs_s} $$

The argument is that the leaves "steal" radiation during the day but somehow work in parallel with the soil during the night. So, if we always use the daytime equation, large leaf area indices will nearly stop evapotranspiration. If we use the second equation instead, large leaf area indices will increase evapotranspiration a little.

During nighttime, condensation often dominates evaporation, and huge leaf area indices enhance condensation. Also, condensed water neither needs to pass the leaves' surface nor the soil's surface. So, wouldn't it be right to always set the surface resistance to zero if other factors suggest the dominance of condensation?

tyralla commented 11 months ago

Surface resistance, part 3

Regarding the last discussion, I think the selected simulation step size should make a difference. With a daily step size (as used by MORECS?), high surface resistances at midnight affect the intense global radiation at noon, whereas an hourly step size (as used by AMBAV?) should prevent such confusion. It is not likely that the same equations fit well in both cases.

tyralla commented 11 months ago

Wind speed / Aerodynamic resistance, part 9

According to Zenker (2003), the equation for aerodynamic resistance I selected above was already used by van Baven (1967), but with a $\kappa = 0.41$:

$$ r_a = \left( \frac{ln(10 / z0)}{0.41} \right)^2 \cdot u{10m} $$

Overall, $\kappa = 0.41$ seems the standard choice in this context, so we should adopt it.

tyralla commented 11 months ago

Surface resistance, part 4

I had a short look at Groh et al. (2019), Monteith (1956), and Iritz & Lindroth (1994). It's all a little contradictory or incomparable (different sites and methods), but all seem to speak for not increasing the combined surface resistance too much during nighttime. So, I choose:

$$ \frac{1}{rs} = \frac{LAI}{2800} + \frac{1}{rs_s} $$

tyralla commented 11 months ago

Surface resistance, part 5

Method Calc_LanduseSurfaceResistance_V1 of Evap-MORSIM includes some special modifications that only address coniferous forest (and not mixed forests). I will not implement such modifications in Evap-AMBAV-1 because AMBAV's documentation does not mention something similar, and I am not aware of a reason to handle coniferous trees so significantly differently from any other plants (which does not mean that there is no good reason).

tyralla commented 11 months ago

Surface resistance, part 6

We discussed the strategy to increase the soil surface resistance when the topmost soil layer dries out. However, we have not found a more convincing way than the simple precipitation index-like approach suggested by the AMBAV documentation. On first thought, one would prefer to directly use something like the relative soil moisture calculated by the main model, but this would be even worse because models like HydPy-L(ARSIM) and HydPy-H(BV) calculate only one average soil moisture value for the whole considered soil column which is a poor predictor for topmost soil layer's wetness.

However, we came up with the idea to not use seasonally varying thresholds that depend on seasonal evaporation patterns but to let the user define a fixed value ($WetnessThreshold$) that must be exceeded by the ratio of the values of precipitation and potential evapotranspiration over the last 24 hours:

$$ r{new} = \begin{cases} r{old} + \Delta &|\ p/e < \tau\ w &|\ p/e \geq \tau \end{cases} $$

$$r = SoilResistance$$

$$w = WetSoilResistance$$

$$ \Delta = SoilResistanceIncrease$$

$$\tau= WetnessThreshold$$

$$p = DailyPrecipitation$$

$$e = DailyPotentialSoilEvapotranspiration$$

This approach complicates the code (we need to remember previous hourly values to calculate the "daily" values) but should significantly improve usability, so it seems worth the effort.

tyralla commented 11 months ago

Albedo, part 2

In our discussions, we decided to introduce four new albedo parameters: SoilAlbedo, SoilAlbedoSnow, LeafAlbedo, and LeafAlbedoSnow. This allows users to choose whether they prefer to handle albedo increases due to snow differently for leaf and soil surfaces (and how).

The mentioned parameters do not allow defining time-variant values (e.g. monthly values), which is consistent with how Evap-AMBAV handles the surface resistance (see above). Hence, there is no need to introduce a new SoilType parameter.

We will reuse the WetnessThreshold approach implemented to adjust the soil's surface resistance.

tyralla commented 11 months ago

Albedo, part 3

I prefer to rename SoilAlbedo and SoilAlbedoSnow to GroundAlbedo and GroundAlbedoSnow because they must also allow defining values for water areas, sealed surfaces, etc.

tyralla commented 11 months ago

Albedo, part 4

Second, it is corrected to 0.8 for complete snow covers and $(\alpha_s + 0.8) / 2$ for incomplete snow covers. (Why does snow not affect the plant albedo?)

Currently, HydPy-L only distinguishes between completely snow-free and snow-covered hydrological response units. However, HydPy-H can provide more detailed information (due to its internal consideration of multiple snow classes per response unit). I propose the following equations to avoid further decisions on what can be considered as "missing", "incomplete", and "complete":

$$ CurrentAlbedo = \alpha_g + (\alpha_l - \alpha_g) \cdot min(lai/4, \ 1) $$

$$ \alpha_g = \omega \cdot GroundAlbedoSnow + (1 - \omega) \cdot GroundAlbedo $$

$$ \alpha_l = \omega \cdot LeafAlbedoSnow + (1 - \omega) \cdot LeafAlbedo $$

$$ \omega = SnowCover $$

$$ lai = LeafAreaIndex $$

tyralla commented 11 months ago

Albedo, part 5

The final version, which also takes the topmost soil layer's wetness into account, is:

$$ CurrentAlbedo = \alpha_g + (\alpha_l - \alpha_g) \cdot min(lai/4, \ 1) $$

$$ \alpha_g = \omega \cdot GroundAlbedoSnow + (1 - \omega) \cdot \alpha_g^* $$

$$ \alpha_g^* = GroundAlbedo \cdot \begin{cases} 1 &|\ \overline{Soil} \ \lor \ p / e < \tau \ 0.5 &|\ Soil \ \land \ p / e \geq \tau \end{cases} $$

$$ \tau = WetnessThreshold $$

$$ p = DailyPrecipitation $$

$$ e = DailyPotentialSoilEvapotranspiration $$

$$ \alpha_l = \omega \cdot LeafAlbedoSnow + (1 - \omega) \cdot LeafAlbedo $$

$$ \omega = SnowCover $$

$$ lai = LeafAreaIndex $$

tyralla commented 11 months ago

Longwave radiation, part 2

We can implement both equations as noted. However, $N$, the sky coverage fraction, and $C$, a coefficient reflecting the cloud type, are often not directly available for simulation.

We discussed using meteorological observations (for calibration) and simulations (for forecasts and predictions) of $N$ and $C$. Reading such time series would be a simple and elegant solution from the model perspective. However, providing this data could cause more effort on the preprocessing side than we can currently spend.

So, we decided to estimate $N$ based on the fraction of the actual and the astronomically possible sunshine duration. During nighttime and for an hourly simulation step size, we will reuse the last simulation step where the possible sunshine duration is one hour but reduce its $N$ value by a factor to account for the often less cloudy conditions during nighttime due to less convection. (On further thinking, the frequent occurrence of fog and low stratus during nights could also indicate setting this factor to values larger than one.)

As we want to keep $N$ up-to-date during the daytime and reuse its last value during the nighttime instead of calculating an average over the previous 24 hours (as we do for other models), adding a state sequence is more straightforward than adding a log sequence. I will name it CloudCoverage.

Regarding $C$, we decided to introduce the plain parameter CloudTypeFactor. We can easily replace it with an equally named input sequence later.

tyralla commented 11 months ago

Longwave radiation, part 3

As we want to keep $N$ up-to-date during the daytime and reuse its last value during the nighttime instead of calculating an average over the previous 24 hours (as we do for other models), adding a state sequence is more straightforward than adding a log sequence. I will name it CloudCoverage.

Things were a little more complicated than expected due to the additional demand to allow for differences between daytime and nighttime cloud cover degrees. I had to split everything into two methods. The first one (Update_CloudCoverage_V1) works as described above:

$$ c_{new} = \begin{cases} s / s_0 &|\ s0 = h \ c{old} &|\ s_0 < h \ \end{cases} $$

$$ c = CloudCoverage $$

$$ s = SunshineDuration $$

$$ s_0 = PossibleSunshineDuration $$

$$ h = Hours $$

The second one (Calc_AdjustedCloudCoverage_V1) performs the desired adjustment and assigns the result to an additional factor sequence:

$$ AdjustedCloudCoverage = \omega \cdot c + (1 - \omega) \cdot n \cdot c $$

$$ \omega = PossibleSunshineDuration / Hours $$

$$ c = CloudCoverage $$

$$ n = NightCloudFactor $$

tyralla commented 11 months ago

Longwave radiation, part 4

$$ RL_{a} = (1 + C \cdot N^2) \cdot \big(1 - 0.261 \cdot exp(-7.77 + 10^{-4} \cdot AirTemperature^2)\big) \cdot \alpha \cdot (AirTemperature + 273.1)^4 $$

There is another error in the given equation. Inserting the albedo neither makes sense nor leads to reasonable values. According to Idso and Jackson (1969) we must instead use the Boltzmann constant:

$$ RL_{a} = (1 + C \cdot N^2) \cdot \big(1 - 0.261 \cdot exp(-7.77 \cdot 10^{-4} \cdot AirTemperature^2)\big) \cdot \sigma \cdot (AirTemperature + 273.1)^4 $$

tyralla commented 11 months ago

Longwave radiation, part 5

In some preliminary experiments, I noticed pronounced differences between the above equation and those of Evap-FAO56 and Evap-MORSIM, which cause reduced evapotranspiration rates. I double-checked its implementation but cannot find an error. We should evaluate this more thoroughly later.

tyralla commented 11 months ago

Saturation vapour pressure

The AMBAV documentation writes the Penman-Monteith equation in two versions—the first deals with specific humidity, and the second deals with vapour pressure. I prefer the second one due to its better fit to the available data and the already implemented methods.

The documentation provides little information regarding the calculation of air humidity-related properties. There is only one equation, and it calculates the saturated specific humidity's slope:

$$ S = \frac{2.54}{(T + 273.3)^2} \cdot exp \left( 19.0785 - \frac{4098.03}{T + 237.3} \right) $$

The temperature delta $273.3$ is also used by the equation underlying method Calc_SaturationVapourPressure_V1 (Allen) and differs from the delta $234.175$ used by method Return_SaturationVapourPressure_V1 (Weischet). So, unless nobody objects, I assume that the methods Calc_SaturationVapourPressure_V1 and Calc_SaturationVapourPressureSlope_V1 are an appropriate choice for Evap-AMBAV.

tyralla commented 11 months ago

Heat capacity of air

The AMBAV documentation states a specific air heat for a constant pressure of 1,005 J/kg/K. We realised that the unit of the corresponding parameter HeatCapacityAir is confusing (WT/kg/K) and agreed on adjusting it to J/kg/K, too.

tyralla commented 11 months ago

Evapotranspiration and snow evaporation

Due to its extraction from the HydPy-L-Land Version 3 and HydPy-L-Land Version 4, HydPy-Evap-MORSIM suppresses evapotranspiration when its main model tells there is a snow cover, which is okay if evaporation during snow periods is negligible or if the snow module calculates evaporation independently from the primary evaporation module (which is the case for Knauf's method).

As discussed above, AMBAV follows its own way for reducing evapotranspiration during snow periods. Hence, I will not include such a suppression mechanism into Evap-AMBAV-1.0 for now. This means, when coupling HydPy-L-Land Version 3 with Evap-AMBAV-1.0, evapo(transpi)ration can be taken into account twice during snow periods.

I suggest to discuss this complex topic separately in another issue and to postpone the eventual adaption of Evap-AMBAV-1.0 until we start with the modularisation of snow processes.

tyralla commented 11 months ago

Preliminary comparison between Evap-MORSIM and Evap-AMBAV-1.0 / Longwave radiation, part 6

I decided to make some preliminary comparisons between both models before finishing Evap-AMBAV-1.0. In the non-tree-vegetation integration test, Evap-MORSIM calculates a potential interception evaporation rate of 4.5 mm/day for a crop height of 10 m. In contrast, Evap-AMBAV-1.0 currently estimates 5.8 mm/day. Why is that?

First finding: The critical difference lies in the aerodynamic resistance. With a reduced crop height of 9.999 m, Evap-MORSIM calculates a rate of 11.1 mm/day (the reason is explained here). So, I guess we can be pretty happy with the current aerodynamic resistance estimates of Evap-AMBAV-1.0 for now, but we should remember this finding if we later conclude that Evap-AMBAV-1.0 models too high evapotranspiration rates for forests.

Second finding: The energy balance calculated by Evap-MORSIM is 138 W/m², while Evap-AMBAV-1.0 calculates a net loss of -41 W/m². Net short-wave radiation is identical in both cases (152 W/m²), and the soil heat flux is relatively similar (3.0 W/m² vs. 5.7 W/m³). So, the critical difference lies in net longwave radiation (23.8 W/m² vs. 192.6 W/m²). This corresponds to the finding reported above. Only Evap-AMBAV-1.0 allows for a clear distinction between earth surface radiation and atmospheric counter radiation, which are 506 W/m² and 314 W/m². The first term is higher than I would expect it. The additional global radiation-related correction term is negligible. So, am I possibly struggling because of another plain typo in the AMBAV documentation? Suppose I replace 293.1 with 273.1, the expected value when applying the Boltzmann law. In that case, outgoing longwave radiation becomes 390 W/m², which reduces net longwave radiation to 76 W/m², which increases total net radiation to 76 (incidentally the same value), which increases potential evapotranspiration to 8.2 mm/day. I think this is the way to go, so if nobody objects, I will replace

$$ RL_{s} = 0.97 \cdot \sigma \cdot (AirTemperature + 293.1)^4 + 0.07 \cdot (1 - \alpha) \cdot GlobalRadiation $$

with

$$ RL_{s} = 0.97 \cdot \sigma \cdot (AirTemperature + 273.1)^4 + 0.07 \cdot (1 - \alpha) \cdot GlobalRadiation $$

tyralla commented 11 months ago

Water areas, part 2

AMBAV does not consider evaporation from water areas. We decided on using the classical Penman equation similar to HydPy-Evap-MORSIM.

First choice(s): We could either (1A) calculate daily evaporation values based on daily input data (averaged over the last 24 hours) or (1B) calculate hourly evaporation values and aggregate them later. 1A agrees with Evap-MORSIM and hence would require little implementation effort. 1A requires less computational effort and would result in a cleaner model (less "log" and "daily" sequences and less methods to calculate them). Additionally, we could decide to enforce a strict averaging to daily values or that users can define a weighting factor to customise the temporal smoothing (like implemented by method Update_EvI_WEvI_V1).

Second choice: We could either (2A) apply the "Dalton wind function" as in method Calc_WaterEvaporation_V3, but then would have to adjust the wind speed to a height of 2 m above the ground or adjust the coefficients. Or (2B), we could use the original Penman equation and insert the aerodynamic resistance for water surfaces.

I choose options 1A (which seems less experimental to me) and 2B (which seems more consistent with the other evapo(transpi)ration calculations performed by Evap-AMBAV-1.0).

tyralla commented 11 months ago

Evaporation from water areas, part 2

I choose options 1A (which seems less experimental to me) and 2B (which seems more consistent with the other evapo(transpi)ration calculations performed by Evap-AMBAV-1.0).

I changed my mind. We would need to calculate additional properties like "daily air density" or "daily aerodynamic resistance". Additionally, it would not always be clear how to calculate them (average the hourly wind speed and directly calculate the daily aerodynamic resistance or average the hourly aerodynamic resistance?). So many considerations and computations for possibly no benefit. So I switch to 1B (with a fixed averaging over the last 24 hours) and 2B.

tyralla commented 10 months ago

Wind speed / Aerodynamic resistance, part 10

In the water-area of Evap-MORSIM, evaporation is 3.2 mm/day. For Evap-AMBAV-1.0, it is 1.9 mm/day.

tyralla commented 10 months ago

Coupling Evap-PET-AMBAV-1.0 to Evap-Minhas

We now have a functioning version of Evap-PET-AMBAV-1.0 with an acceptable set of integration tests but no main model we can couple it with, which hinders us from performing the first real-work test runs. As it is still unclear if we want to implement Evap-AET-AMBAV-1.0 now or later, we decided on making Evap-Minhas and Evap-AET-HBV96 ready for this task. This requires additional thought because Evap-PET-AMBAV-1.0 and the already implemented PET models follow different interfaces. The "old" interface assumes only one kind of potential evapotranspiration, and hence (in short version) looks so:

class PETModel_V1(modeltools.SubmodelInterface):

    typeid: ClassVar[Literal[1]] = 1

    def determine_potentialevapotranspiration(self) -> None:
    def get_potentialevapotranspiration(self, k: int) -> float:
    def get_meanpotentialevapotranspiration(self) -> float:

Whereas the new interface for Evap-AMBAV-1.0 looks so:

class PETModel_V2(modeltools.SubmodelInterface):

    typeid: ClassVar[Literal[2]] = 2

    def determine_potentialinterceptionevaporation(self) -> None:
    def get_potentialinterceptionevaporation(self, k: int) -> float:
    def determine_potentialsoilevapotranspiration(self) -> None:
    def get_potentialsoilevapotranspiration(self, k: int) -> float:
    def determine_potentialwaterevaporation(self) -> None:
    def get_potentialwaterevaporation(self, k: int) -> float:

This difference is not only of technical but also hydrological importance because neither Evap-Minhas nor Evap-AET-HBV96 have any means to work with different types of evapotranspiration so far.

My first idea regarding potential interception evaporation and soil evapotranspiration is to split the old "combined" potential interception evapotranspiration sequence into two. Then, we had two tasks:

  1. Ensure that Evap-Minhas and Evap-AET-HBV96 query the same value when working with PETModel_V1 but different values when working with PETModel_V2 (the technical part)
  2. Ensure that the process equations of Evap-Minhas and Evap-AET-HBV96 work in cases where the given potential interception evaporation and soil evapotranspiration values differ (the hydrological part)

This requires a little more work and is more disruptive than expected. Still, it separates everything cleanly and helps track the calculation of the different PET components in later applications.

tyralla commented 10 months ago

Coupling Evap-PET-AMBAV-1.0 to Evap-Minhas, part 2

This is how Evap-Minhas currently calculates soil evapotranspiration:

$$ SoilEvapotranspiration = (pet - e_i) \cdot \frac{1 - exp\left(-f\cdot \frac{w}{m } \right)}{1 + exp\left(-f\cdot \frac{w}{m } \right) - 2 \cdot exp(-f)} $$

$$ pet= PotentialEvapotranspiration $$

$$ e_i = InterceptionEvaporation $$

$$ f = DisseFactor $$

$$ w = SoilWater $$

$$ m = MaxSoilWater $$

The mathematically obvious generalisation for different values of potential interception evaporation ($pe_i$) and potential soil evapotranspiration is ($pet_s$):

$$ SoilEvapotranspiration = pet_s \cdot \left( 1 - \frac{e_i }{pe_i} \right) \cdot \frac{1 - exp\left(-f\cdot \frac{w}{m } \right)}{1 + exp\left(-f\cdot \frac{w}{m } \right) - 2 \cdot exp(-f)} $$

I will adjust Evap-Minhas accordingly if no one suggests a better equation.

I just saw that the changed term corresponds to Update_SoilEvapotranspiration_V3 (the "Wigmosta" approach). I will reuse this method to avoid code duplications.

tyralla commented 10 months ago

Coupling Evap-PET-AMBAV-1.0 to Evap-AET-HBV96

Evap-AET-HBV96 already implements three methods for calculating and updating soil evapotranspiration.

Method Calc_SoilEvapotranspiration_V1 applies the standard HBV "LP" relationship:

$$ et_s = pet \cdot \frac{s}{l \cdot m} $$

$$ et_s = SoilEvapotranspiration $$

$$ pet = PotentialEvapotranspiration $$

$$ s = SoilWater $$

$$ l = SoilMoistureLimit $$

$$ m = MaxSoilWater $$

Next, method Update_SoilEvapotranspiration_V1 performs an excess reduction (meaning, it is only applied if there is an excess), which is the current critical point for us:

$$ et_s^{new} = et_s^{old} - r \cdot (et_s^{old} + e_i - pet) $$

$$ r = ExcessReduction $$

$$ e_i = InterceptionEvaporation $$

Finally, method Update_SoilEvapotranspiration_V2 reduces soil evapotranspiration due to possible snow covering, which can stay as is:

$$ et_s^{new} = (1 - c) \cdot et_s^{old} $$

$$ c = SnowCover $$

The first adjustment is trivial:

$$ et_s = pet_s \cdot \frac{s}{l \cdot m} $$

$$ pet_s = PotentialSoilEvapotranspiration $$

The intention of the first update method is to restrict $et_s + e_i$ to $pet$ (for $r = 1$), to $2 \cdot pet$ (for $r = 0$), or something in between. My suggestion for $r = 0$ is to restrict $et_s + e_i$ to $pet_s + pe_i$ and for $r = 1$ to $pe_i$ (the latter because interception evaporation usually takes priority). So, the basic equation could be:

$$ et_s^{new} = et_s^{old} - r \cdot \Big(et_s^{old} + e_i - \big(r \cdot pe_i + (1 - r) \cdot (pet_s + pe_i) / 2 \big) \Big) $$

One can verify that this equation does not change the HBV96-like behaviour by setting $pet_s = pe_i = pet$ and simplifying.

There might be some additional corner cases (e.g. different signs of $pet_s$ and $pe_i$) which might need special attention. I will try to cover them all in the unit tests and only report them here if their solution is not apparent.

tyralla commented 6 months ago

It's time to complete HydPy-PET-AMBAV-1.0 and close this issue. There are multiple uncertainties, but maybe we can narrow them down to the most important ones:

Currently, HydPy-PET-AMBAV-1.0 calculates surprisingly low evaporation values for water areas and high evapotranspiration values from tree-like vegetation. One gets a good impression of this by comparing the integration tests of evap_morsim and evap_pet_ambav1.

Regarding evaporation from water areas, the relevant difference seems to be that MORSIM uses a version of the Dalton wind function for calculating surface resistance (DVWK 1996), while AMBAV uses the original Penman-Montheith equation with zero surface resistance. We could either use the Dalton approach also for AMBAV or just modify the roughness length, which is calculated based on the crop height. I performed a sensitivity analysis and achieved highly similar results in the mentioned integration tests by setting the "crop height" of water areas to 1 m.

Regarding evapotranspiration from tree-like vegetation, I suppose the underlying equations for calculating the aerodynamic resistance of both models do not work well for forests. However, MORSIM trims the resistance of crops taller than 10 m to a specific value, while AMBAV does not. I don't know if the fixed MORSIM value is well justified. If so, we could use it in AMBAV, as well.

tyralla commented 6 months ago
  • Change the calculation of evaporation for water areas?

I played around and maybe I was wrong. When letting HydPy-PET-AMBAV-1.0 use the Dalton approach, the evaporation from the water area increased to a lesser extent than expected. The difference in the net longwave radiation estimates seems to be the major cause of the evaporation differences. I now tend to say we require deeper investigations and see how AMBAV works in real applications before doing any ad-hoc changes.

tyralla commented 6 months ago
  • Change the calculation of the aerodynamic resistance for tree-like vegetation?

I did a web search but could not find much information on estimating aerodynamic resistance for forests (based only on wind speed measurements at a weather station somewhere else). I suggest reusing the approach already implemented in HydPy-MORSIM but with a slight modification.

https://github.com/hydpy-dev/hydpy/issues/111#issuecomment-1515826753 shows that aerodynamic resistance decreases with increasing crop height until 10 m is reached. Then, aerodynamics suddenly rises and stays at the same level for all higher crop heights.

Maybe we should at least ensure a continuous transition between both equations. According to our documentation, we could turn the switch at a crop height of 1.14 m instead of 10 m. We can accomplish this by extending the calculation of the AerodynamicResistanceFactor:

$$ max \big( ln(2 / z_0) \cdot ln(10 / z_0) / 0.41^2, \ 94 \big) $$

The hard-coded value "94" might be too restrictive. So, I will introduce a "fixed parameter" (AerodynamicResistanceFactorMinimum), which allows us to experiment with alternative values later.

tyralla commented 6 months ago

I see that AMBAV suggests reducing the wind speed at a height of 10 m by a factor of 0.6 above trees. We already implemented this. So, following my last idea, we would solve the same problem twice and likely "overcorrect" the overestimation of evapotranspiration from trees.

We could remove the wind speed adjustment or set the threshold parameter to $0.6 \cdot 94$. The first means a little more effort but is cleaner, so I prefer it.

tyralla commented 6 months ago

I changed HydPy-PET-AMBAV-1.0 as explained. I now consider it operational and merge it into the master branch after the last tests. If we find the need to improve some nuances later, we should open a specific new issue for that.