Support using GARTO as a submodel of HydPy-L

[tyralla]

We currently aim to support LARSIM-like simulations with HydPy-L that better take infiltration excess into account via a Green & Ampt (GA) infiltration approach. However, HydPy-L is already very complex, and a GA approach suitable for continuous simulations is not too simple either. So, combining the equations of both models is something we should avoid. Instead, we strive for an independently usable GA model that allows for coupling with HydPy-L (and, at best, other "complete" hydrological models like HydPy-H). Unfortunately, this is beyond HydPy's current modularisation functionalities. Extending these functionalities will be non-trivial, but we think it will pay off very soon due to ways to keep our application models smaller and increase the users' flexibility in combining different model modules.

The necessary work will touch on many topics and require many decisions. Hence, we split the discussion into three issues that separately deal with the GA methodology itself (#89), the general approach to couple "complete" application models like HydPy-L with "submodels" like GA (#90), and the specific coupling interface for HydPy-L and GA (#91).

[tyralla]

For now, I will name the required interface SoilInterface. (First, I thought about InfiltrationInterface, but the main model and its submodel must also exchange evaporation data (and, eventually, percolation data).

[tyralla]

For coupling HydPy-L with HydPy-GA, I would think of the following workflow (for a start, I neglect capillary rise and possibly some other special cases):

1. L calculates WaDa (water reaching the soil routine) as usual.
2. L sends WaDa to GA as InitialSurfaceWater.
3. L calculates QDB (direct runoff release from the soil storage (misleading description)) as usual.
4. L sends WaDa - QDB to GA as ActualSurfaceWater.
5. L calculates EvB (actual evaporation of soil water) as usual.
6. L sends EvB to GA as Demand.
7. L triggers GA's Perform_GARTO_V1 method.
8. L queries GA's Infiltration to increase QDB eventually.
9. L queries GA's Withdrawal and takes it as the new EvB.
10. L queries GA's Percolation and takes it as the new QBB.
11. L calculates QIB1 and QIB2 (the first and the second component of the interflow release from the soil storage) as usual.
12. L sends QIB1 + QIB2 to GA as Demand.
13. L triggers GA's Withdraw_AllBins_V1 method.
14. L updates its soil water content (preferably by querying it from GA?).

[tyralla]

Things I do not like or feel unsure about the above workflow:

• Evaporation withdrawal occurs during GA's numerical substeps (when calling Perform_GARTO_V1), but interflow withdrawal occurs at the end of the simulation interval.
• Splitting L's surface water supply into initial and actual surface water depth seems smart at the moment, but we still need to check if it works.
• Is GA's percolation a good surrogate for L's groundwater recharge generation?
• L calculates QDB, QIB1, QIB2, QBB, EvB, and QKap in parallel (based on the same initial soil water content) and uses them to update the soil water content in one step. Very early during the implementation of L, we thought this was what LARSIM does (or did), but we are not so sure anymore. Howsoever, this state updating after parallel flux calculations can complicate things (especially when we need to respect thresholds like zero water content). I have a gut feeling that it will complicate coupling with GA, too.

[tyralla]

• Evaporation withdrawal occurs during GA's numerical substeps (when calling Perform_GARTO_V1), but interflow withdrawal occurs at the end of the simulation interval.

Possible solutions:

1. Withdraw evaporation also just once, at the end of the simulation interval.
2. Let GA calculate both interflow components on its own.

The firsts solution could go for a workflow like this:

1. L calculates WaDa (water reaching the soil routine) as usual.
2. L sends WaDa to GA as InitialSurfaceWater.
3. L calculates QDB (direct runoff release from the soil storage (misleading description)) as usual.
4. L sends WaDa - QDB to GA as ActualSurfaceWater.
5. L triggers GA's Perform_GARTO_V1 method.
6. L queries GA's Infiltration to increase QDB eventually.
7. L queries GA's Percolation and takes it as the new QBB.
8. L calculates EvB (actual evaporation of soil water) as usual.
9. L calculates QIB1 and QIB2 (the first and the second component of the interflow release from the soil storage) as usual.
10. L sends EvB + QIB1 + QIB2 to GA as Demand.
11. L triggers GA's Withdraw_AllBins_V1 method.
12. L queries GA's Withdrawal and uses it to reduce EvB, QIB1, or QIB2 when necessary.
13. L updates its soil water content (preferably by querying it from GA?).

Before thinking about the second point, we maybe need a clearer picture of the responsibilities of both models.

[tyralla]

Regarding responsibility: So far, I have thought about handling the water content in both models (BoWa in L and Moisture times FrontDepth in GA). But it's preferable if only one model needs to handle it, and it must be GA because it relies on more detailed information. This means we could remove the "state sequence" BoWa (or turn it into an "aide sequence" for convenience) and let L use something like the current watercontent property of GA instead.

I try to combine this idea with the above workflow in the form of pseudo-code:

1. bowa = ga_model.calc_watercontent()
2. wada = l_model.calc_wada()
3. ga_model.set_initialsurfacewater(wada)
4. qdb = l_model.calc_qdb(bowa)
5. ga_model.set_actualsurfacewater(wada - qdb)
6. ga_model.perform_garto()
7. qdb += (wada - qdb) - ga_model.get_infiltration()
8. qbb = ga_model.get_percolation()
9. # bowa = ga_model.calc_watercontent()  update bowa? maybe before each usage?
10. evb = l_model.calc_evb(bowa)
11. qib1 = l_model.calc_qib1(bowa)
12. qib2 = l_model.calc_qib2(bowa)
13. ga_model.set_demand(evb + qib1 + qib2)
14. ga_model.withdraw_allbins()
15. (evb, qib1, qib2) *= ga_model.get_withdrawal() / (evb + qib1 + qib2)

It appears a little random that GA calculates partly QDB and (partially?) QBB but not at all QIB1 and QIB2. However, moving everything to GA would eventually make coupling harder (for example, HBV has different ideas about interflow than LARSIM). So, I think it is okay to start like this. Maybe we can make further modularisation improvements later by making the interflow calculation a submodel of GA?

[tyralla]

On the capillary rise issue.

Just for completeness: The original GARTO model considers dynamical groundwater fronts. But going this way for replacing LARSIM's capillary rise would be much too much effort and increase computation times severely.

So, if we calculate the capillary rise in the LARSIM way, how (or where) do we add it to GA's water content? Ideas:

1. Let a groundwater front rise. (already rejected above)
2. Add an "inverse method" of Withdraw_AllBins_V1, something like Supply_AllBins_V1, which goes from the left to the right bin (from low to high moisture) instead from right to left. (Then, we might need to change the current meaning of "supply" related to the surface water depth - but it seems we need to change the related method anyhow.)
3. Let Withdraw_AllBins_V1 support negative withdrawals? (This idea does not seem to make much sense.)

So, the second idea seems to be the favourite. The pseudo-code then could look like this:

1. bowa = ga_model.calc_watercontent()
2. wada = l_model.calc_wada()
3. ga_model.set_initialsurfacewater(wada)
4. qdb = l_model.calc_qdb(bowa)
5. ga_model.set_actualsurfacewater(wada - qdb)
6. ga_model.perform_garto()
7. qdb += (wada - qdb) - ga_model.get_infiltration()
8. qbb = ga_model.get_percolation()
9. # bowa = ga_model.calc_watercontent()  update bowa? maybe before each usage?
10. qkap = l_model.calc_qkap(bowa, groundwater)
11. ga_model.set_supply(qkap)
12. ga_model.supply_allbins()
13. qkap = ga_model.get_addition()  # for the unlikely case capillary rise results in a soil moisture overflow
14. evb = l_model.calc_evb(bowa)
15. qib1 = l_model.calc_qib1(bowa)
16. qib2 = l_model.calc_qib2(bowa)
17. ga_model.set_demand(evb + qib1 + qib2)
18. ga_model.withdraw_allbins()
19. (evb, qib1, qib2) *= ga_model.get_withdrawal() / (evb + qib1 + qib2)

[tyralla]

Under the above discussion, a prototype of our SoilInterface (or should we name it SoilRoutine, which sounds more hydrological than informatical?):

class SoilInterface:
    def set_initial_surface_water(...): ...
    def set_actual_surface_water(...): ...
    def set_soil_water_supply(...): ...
    def set_soil_water_demand(...): ...
    def execute_infiltration(...): ...
    def add_soil_water(...): ...
    def remove_soil_water(...): ...
    def get_infiltration(...): ...
    def get_percolation(...): ...
    def get_soil_water_addition(...): ...
    def get_soil_water_removal(...): ...
    def get_soil_water_content(...): ...

A technical issue: We can generally assume to deal with a vector of soil columns. If the main model handles a single soil (like WALRUS) or a 2-dimensional soil grid (not implemented so far), it can do the necessary conversion (e.g. flattening) by itself. But we must decide if we want to pass complete vectors or single vector entries. I think the letter often comes with more boilerplate code. However, it might often be more straightforward, especially concerning avoiding computational overhead due to array creations or shape modifications. So, I would start with passing individual vector entries.

A conceptual question: How many main models can use such an interface, and how many submodels fit into it?

[tyralla]

A very preliminary (pure Python) prototype of our coupling of HydPy-L and HydPy-GA is running. Before continuing technically, the numerous decisions required for making the coupling hydrologically sound need attention, especially regarding the different used volume-related soil parameters. As things are pretty complicated, we first restrict the discussion to lland_v2 and ga_garto_sub1 (should also apply to lland_v1, but lland_v3 and lland_v4 need special consideration due to using the Penman-Monteith equation) and neglect capillary rise.

ga_garto_sub1 defines the following volume-related soil parameters:

• ResidualMoisture $\theta_r$
  • There is no (liquid) water flow below this value (actual conductivity is zero).
  • There is no evaporation below this value.
  • Hence, the actual moisture becomes never smaller than $\theta_r$.
• SaturationMoisture $\theta_s$
  • At this point, actual conductivity reaches saturated conductivity.
  • At this point (given enough rainfall), infiltration, percolation, and saturated conductivity are identical.
  • Hence, the actual moisture becomes never larger than $\theta_s$.

A concept like "field capacity". So, besides minor numerical inaccuracies, there is no generation of percolation during zero-rainfall periods.

lland_v2 defines the following volume-related soil parameters:

• PWP (permanent wilting point)
  • Does not affect evapotranspiration despite its name.
  • If RBeta is false, PWP serves as the starting point for percolation (groundwater recharge)
  • If RBeta is true, there is no percolation below FK (which prevents groundwater recharge and capillary rise from coinciding).
  • Below PWP, the generation of the slow interflow component is always zero.
• FK (field capacity)
  • If RBeta is true, FK is the starting point for percolation
  • If RBetta is false, below FK, the relation between soil moisture and percolation is linear.
  • Either way, above FK, there can be a nonlinear increase of percolation with soil moisture.
  • FK is the starting point for the generation of the fast interflow component.
• WMax (maximum water storage)
  • The highest possible value of soil moisture.
  • If soil moisture reaches WMax, actual soil evapotranspiration reaches its highest possible value (potential soil evaporation).
  • If soil moisture reaches WMax, actual percolation reaches its highest possible value ( $Beta \cdot FBeta $).
  • If soil moisture reaches WMax, direct runoff generation reaches its highest possible value (WaDa).
  • If soil moisture reaches WMax, slow interflow generation reaches its highest possible value (DMin).
  • If soil moisture reaches WMax, fast interflow generation reaches its highest possible value ( $Dmax - DMin$).

possible configuration

1. Get SoilDepth, SaturationMoisture, FK, and ResidualMoisture from somewhere.
2. Assume $SaturationMoisture \cdot SoilDepth = WMax$ (Would it be helpful if the interface calculates WMax automatically? Hard to generalise, but perhaps not impossible.)
3. Set DMin to zero (now standard practice in LARSIM and simplifies understanding PWP)
4. Set PWP to FK (so that no percolation occurs below FK).
5. Set RBeta to true (just for documentation; the real value does not really matter when setting PWP to FK).
6. Set Beta and FBeta to sufficiently high values so that an oversaturated soil drains within two or three days to FC. (Are there approaches to relating this drainage's speed to soil properties? It would be great to avoid calibrating Beta and FBeta in the future).
7. Use DMax for calibrating the ratio between interflow and baseflow.

interface details and effects

Going this way, it can happen for arid conditions that lland_v2 calculates small amounts of actual soil evaporation, but ga_garto_sub1 cannot withdraw this amount from the soil. However, I think this slight inconsistency is negligible, but we need to correct the actual soil evaporation when necessary so that the water balance stays intact.

As discussed elsewhere, lland_v2 withdraws interflow from ga_garto_sub1 as it withdraws evaporation.

lland_v2 would calculate percolation due to oversaturation of the soil and ga_garto_sub1 as a direct response to rainfall. Both functionalities are reasonable so I would take the sum of both percolation components as groundwater recharge.

[tyralla]

4. qdb = l_model.calc_qdb(bowa)
5. ga_model.set_actualsurfacewater(wada - qdb)
6. ga_model.perform_garto()
7. qdb += (wada - qdb) - ga_model.get_infiltration()

If we want to combine percolation from land_v2 and ga_garto_sub1, we must withdraw the percolation component calculated by lland_v2 from the soil water content of ga_garto_sub1.

[tyralla]

A preliminary comparison.

The results of our first integration test of lland_v1:

This is what it looks like when coupling lland_v1 with ga_garto_sub1, assuming the same loam soil as we are using in the integration tests of ga_garto:

Clearly, rainfall intensity now plays a much more significant role. I could not detect any water balance errors so far (unless activating capillary rise, which we still need to support).

Next steps (regarding this issue):

• Build a Cython prototype
• Support the capillary rise option
• Solve the name conflict for "supply"
• Finalise the stand-alone model (ga_garto)
• Finalise the submodel (ga_garto_sub1)
• Finalise the main models (lland_v1, lland_v2, lland_v3, lland_v4)

[tyralla]

To test if lland_v3 uses the ga_garto_submodel1 correctly, I duplicated the following test, where a constant rainfall intensity of 20 mm/d wettens up a storage volume of 309 mm over less than a month:

After coupling with ga_garto_submodel1, soil moisture does not exceed 220 mm anymore:

After a few thoughts, I came to the following explanation:

When comparing the |lland_v1| examples :ref:`lland_v1_acker_summer` and
:ref:`lland_v1_acker_garto`, we see an increase in direct runoff generation due to
including infiltration excess when using |ga_garto_submodel1|.  However, comparing the
following results with the :ref:`lland_v3_acker_heavy_rain_daily` example, we see a
decrease in direct runoff generation. This behaviour results from the lower rainfall
rates, which never exceed the soil's (saturated) conductivity. Hence,
|ga_garto_submodel1| always lets the surface water supplied by |lland_v3| infiltrate
and creates no additional direct runoff component. But, starting from August 13, it
calculates extra groundwater recharge. At this time, the (partly saturated) wetting
front reaches the soil's bottom and enables the continuous percolation of all rainfall
through the soil column:

So I think everything is good, and it's nice to have this example in the test suite and documentation. Still, this little surprise shows that we must be cautious when applying the coupling. And we should not forget to consider alternative lower boundary conditions.