HydPy-GA-GARTO (main model for calculating the Green-Ampt / Talbot-Ogden infiltration with redistribution)

ga_garto implements GARTO, a “Green-Ampt infiltration with Redistribution” model that “incorporates features from the Talbot-Ogden infiltration and redistribution method”, developed by Lai et al. (2015). Thankfully, two of its authors, Fred L. Ogden and Cary A. Talbot, provided their source code (written in C), which we used to ensure ga_garto mostly works like this original implementation. However, there are also some differences, which interested users find discussed in issue 89. The most notable difference is that ga_garto generally assumes a hydrostatic groundwater table, while the complete GARTO method also simulates groundwater dynamics (in terms of vertically moving groundwater fronts).

GARTO combines the oriinal Green-Ampt infiltration model with redistribution equations. In the Green-Ampt model, a saturated wetting front moves downwards in a piston-like form through a homogeneous soil with constant initial moisture and assumes rainfall intensity always exceeds infiltration intensity. The redistribution equations apply for periods with little or no rainfall. When redistribution begins, a previously saturated wetting front still grows in length but reduces its relative water content. With the next heavy rainfall event, a new saturated front develops, and both the old and the new front advance downwards.

One improvement of GARTO over simpler GAR approaches is that the number of active wetting fronts is (principally) not restricted. Wetting fronts only disappear when being merged with other fronts, which happens when a newer (wetter) front overshoots an older (dryer) one or when a wetting front reaches the soil’s depth (the groundwater table). Due to approximating natural soil moisture distributions with multiple wetting fronts, GARTO can cope with complex rainfall patterns with (long) hiatus periods, making it a suitable tool for long-term simulations.

Another functionality required for long-term simulations is evapotranspiration, which is not discussed by Lai et al. (2015) but implemented in the mentioned C code. ga_garto requires the complete withdrawal to be calculated externally and takes it from the surface water and the wettest wetting fronts like one can do with the author’s GARTO code for evaporation. The additional transpiration features of the C code (e.g. taking rooted depth into account) are neglected since we expect ga_garto to be applied for modelling the rooted soil zone.

ToDo: Taking the rooted depth as the depth of the considered soil domain and assuming a

hydrostatic groundwater table at the soil’s bottom do not go hand in hand for many catchments. Do we need alternative lower boundary conditions?

Integration tests

Note

When new to HydPy, consider reading section Integration Tests first.

We prepare an initialisation period of 24 hours but restrict the actual simulation and evaluation period to 5 hours for the first examples:

>>> from hydpy import pub
>>> pub.timegrids = "2000-01-01 00:00", "2000-01-02 00:00", "1h"
>>> pub.timegrids.sim.lastdate = pub.timegrids.eval_.lastdate = "2000-01-01 05:00"

We prepare a ga_garto model, set the base time for defining (most) time-dependent parameter values to one hour, connect the model to an Element instance, and hand this element to an IntegrationTest instance that will execute and evaluate the following example simulations:

>>> from hydpy.models.ga_garto import *
>>> parameterstep("1h")
>>> from hydpy import Element, IntegrationTest
>>> soil = Element("soil")
>>> soil.model = model
>>> test = IntegrationTest(soil)
>>> test.dateformat = "%H:%M"
>>> test.plotting_options.axis1 = (inputs.rainfall, fluxes.infiltration,
...                                fluxes.surfacerunoff, fluxes.percolation,
...                                fluxes.withdrawal)
>>> test.plotting_options.axis2 = states.frontdepth

For the first set of examples, we initialise a single soil compartment (via parameter NmbSoils) consisting of three “bins” (via parameter NmbBins). So, besides the first “filled” bin, there can be at most two active wetting fronts at the same time:

>>> nmbsoils(1)
>>> nmbbins(3)

We set the length of the numerical substeps to one second:

>>> with pub.options.parameterstep("1s"):
...     dt(1.0)

The single soil compartment is not sealed, meaning infiltration is possible:

>>> sealed(False)

The defined subarea is irrelevant as long as we deal with a single compartment:

>>> soilarea(2.0)

We perform all tests for an “average” soil (loam) and take its parameters from Lai et al. (2015), table 1:

>>> residualmoisture(0.027)
>>> saturationmoisture(0.434)
>>> saturatedconductivity(13.2)
>>> poresizedistribution(0.252)
>>> airentrypotential(111.5)

In agreement with the evaluations performed by Lai et al. (2015), all our test runs start with a relative moisture content of 10 %:

>>> test.inits = ((states.moisture, 0.1),
...               (states.frontdepth, 0.0),
...               (logs.moisturechange, 0.0))

5 hours, no evaporation, no capillary rise

In this subsection, we focus on a short simulation period of five hours without any evaporation losses or capillary rise gains. The chosen rainfall pattern stems from Lai et al. (2015), table 2 (soil number 4, loam):

>>> inputs.rainfall.series[:5] = [40.0, 0.0, 0.0, 40.0, 0.0]
>>> inputs.evaporation.series = 0.0
>>> inputs.capillaryrise.series = 0.0

The next commands store the current initial conditions, so we can use them later to check that ga_garto is not violating the water balance:

>>> test.reset_inits()
>>> conditions = model.conditions

deep soil

In this example, we set the soil compartment’s depth to one meter so that it does not restrict the movement of any front:

>>> soildepth(1000.0)

The infiltration sums for the first and the second event agree well but not perfectly with the results reported by Lai et al. (2015), table 3 (soil number 4, loam), which are 38.62 mm (instead of 38.90 mm) and 29.67 mm (instead of 30.01 mm). We checked if the different lengths of the numerical substeps (10 s vs 1 s) cause this error, but they do not. The authors’ GARTO code, which we adapted to the same example, also using a numerical substep length of 1 s, calculates the same infiltration sums as ga_garto (38.90 mm and 30.01 mm).

The water content equals the saturation water content at the end of both rainfall impulses and decreases afterwards. Inspecting the front depths’, it looks like everything takes place in the second bin. However, ga_garto actually creates a second wetting front in the third bin during the second impulse, but this front reaches the first one during the same simulation step, so they both get merged in the second bin:

>>> test("ga_garto_5h_1000mm")
Click to see the table
Click to see the graph

There is no indication of an error in the water balance:

>>> from hydpy import round_
>>> round_(model.check_waterbalance(conditions))
0.0

200 mm depth

Next, we repeat the above example with a reduced soil depth of 200 mm:

>>> soildepth(200.0)

All results for the first three hours are identical to those of the deep soil example. But during the second rainfall event, the saturated wetting front reaches the soil’s bottom so that the first “filled” bin becomes saturated and is the only active bin from then on. Complete saturation comes with a dramatic increase in actual conductivity through the whole soil column and thus in percolation (groundwater recharge). On the other hand, the overall downward water movement and thus infiltration slows due to losing capillary drive after building complete contact with groundwater:

>>> test("ga_garto_5h_200mm")
Click to see the table
Click to see the graph

There is no indication of an error in the water balance:

>>> round_(model.check_waterbalance(conditions))
0.0

150 mm depth

When reducing the soil depth further to 150 mm, the first (and, so far, only) wetting front reaches the soil’s bottom yet in the first redistribution phase within the third hour. The soil water deficit at the start of the second rainfall event is smaller due to increased relative moisture and the reduced soil depth, so infiltration declines even more than in the 200 mm depth example:

>>> soildepth(150.0)
>>> test("ga_garto_5h_150mm")
Click to see the table
Click to see the graph

There is no indication of an error in the water balance:

>>> round_(model.check_waterbalance(conditions))
0.0

100 mm depth

When reducing the soil depth even to 100 mm, the soil already saturates in the first hour. Then, the (potential) infiltration sum of the second event gets its smallest possible value defined by saturated conductivity:

>>> soildepth(100.0)
>>> test("ga_garto_5h_100mm")
Click to see the table
Click to see the graph

There is no indication of an error in the water balance:

>>> round_(model.check_waterbalance(conditions))
0.0

5 hours, evaporation

In this subsection, we repeat some of the above examples with an evaporation rate of 1 mm/h:

>>> inputs.evaporation.series = 1.0

deep soil

Repeating the deep soil example but including evaporation reduces surface runoff generation. This effect emerges especially for the second event, where evaporation has a cumulative impact due to also withdrawing water during the previous redistribution phase:

>>> soildepth(1000.0)
>>> test("ga_garto_5h_1000mm_evap")
Click to see the table
Click to see the graph

There is no indication of an error in the water balance:

>>> round_(model.check_waterbalance(conditions))
0.0

100 mm depth

When repeating the 100 mm depth example with evaporation, ga_garto does not always withdraw water from the surface or the currently wettest wetting front (as in the deep soil example) but (during rain-free periods) also from the filled bin:

>>> soildepth(100.0)
>>> test("ga_garto_5h_100mm_evap")
Click to see the table
Click to see the graph

There is no indication of an error in the water balance:

>>> round_(model.check_waterbalance(conditions))
0.0

4 mm depth

Reducing the soil depth to 4 mm is surely no realistic configuration, but we use it to show ga_garto works stably in case the soil runs completely dry during the third hour:

>>> soildepth(4.0)
>>> test("ga_garto_5h_4mm_evap")
Click to see the table
Click to see the graph

There is no indication of an error in the water balance:

>>> round_(model.check_waterbalance(conditions))
0.0

deep soil, 10 sec

In the following, we reuse the deep soil example to give an impression of the relation between computational effort and numerical accuracy. Please consider this the first impression because such relationships can be distinctly different for other conditions. For example, setting the numerical substep length to 10 s (as done by Lai et al. (2015)) instead of 1 s (as in all our previous examples) increases the total infiltration sum from 69.71 mm to 69.75 mm, which is an irrelevant deviation for practical applications:

>>> soildepth(1000.0)
>>> with pub.options.parameterstep("1s"):
...     dt(10)
>>> derived.nmbsubsteps.update()
>>> test("ga_garto_5h_1000mm_evap_10sec")
Click to see the table
Click to see the graph

There is no indication of an error in the water balance:

>>> round_(model.check_waterbalance(conditions))
0.0

deep soil, 1 min

Lai et al. (2015) recommend using a numerical substep length of less than 1 min. If we set it to 1 min exactly, the total infiltration sum increases from 69.7 mm to 69.9 mm. Hence, at least for deep loam soils, calculating in 1 min intervals seems acceptable for practical applications:

>>> with pub.options.parameterstep("1m"):
...     dt(1)
>>> derived.nmbsubsteps.update()
>>> test("ga_garto_5h_1000mm_evap_1min")
Click to see the table
Click to see the graph

There is no indication of an error in the water balance:

>>> round_(model.check_waterbalance(conditions))
0.0

deep soil, 1 hour

This example shows that even when selecting way too long substeps (in this case, equal to the simulation step size), ga_garto works stably and holds the water balance (still, its results might be more or less random):

>>> with pub.options.parameterstep("1h"):
...     dt(1.0)
>>> derived.nmbsubsteps.update()
>>> test("ga_garto_5h_1000mm_evap_1h")
Click to see the table
Click to see the graph

There is no indication of an error in the water balance:

>>> round_(model.check_waterbalance(conditions))
0.0

5 hours, capillary rise

In this subsection, we repeat some of the above examples with a capillary rise rate of 1 mm/h:

>>> inputs.evaporation.series = 0.0
>>> inputs.capillaryrise.series = 1.0

We set the numerical substep length to 1 s, as in most examples above:

>>> with pub.options.parameterstep("1s"):
...     dt(1)
>>> derived.nmbsubsteps.update()

deep soil

Repeating the deep soil example but including capillary rise reduces surface runoff generation. This effect emerges especially for the second event, where capillary rise has a cumulative impact due to also adding water during the previous redistribution phase:

>>> soildepth(1000.0)
>>> test("ga_garto_5h_1000mm_caprise")
Click to see the table
Click to see the graph

There is no indication of an error in the water balance:

>>> round_(model.check_waterbalance(conditions))
0.0

100 mm depth

After decreasing the soil’s depth to 100 mm, the soil can absorb the available capillary rise only at the start of the simulation. The actual soil water addition is 0.82 mm for the first hour instead of the potential 1.0 mm. Afterwards, the soil is full and cannot absorb any capillary rise:

>>> soildepth(100.0)
>>> test("ga_garto_5h_100mm_caprise")
Click to see the table
Click to see the graph

There is no indication of an error in the water balance:

>>> round_(model.check_waterbalance(conditions))
0.0

capillary rise and evaporation

In this example, we show what happens when simultaneously considering both capillary rise and evaporation by assuming the same intensity of 1 mm/h:

>>> inputs.evaporation.series = 1.0

At first, one might assume ga_garto calculates the same infiltration rates as in the deep soil example because the capillary rise and evaporation cancel each other out. However, we see slightly higher infiltration rates as in the deep soil results due to the different “priorities” of both processes. For evaporation, method Withdraw_AllBins_V1 prefers to take water from the wettest bin (proceeds from right to left). For capillary rise, method Withdraw_AllBins_V1 prefers adding water to the dryest bin (proceeds from left to right):

>>> soildepth(1000.0)
>>> test("ga_garto_5h_1000mm_evap_caprise")
Click to see the table
Click to see the graph

There is no indication of an error in the water balance:

>>> round_(model.check_waterbalance(conditions))
0.0

24 hours, no evaporation

In this and the following subsection, we deal with more complex series of rainfall events spanning 24 hours:

>>> pub.timegrids = pub.timegrids.init

The examples demonstrate how ga_garto works when pronounced hiatus periods occur, where the accuracy of the redistribution process becomes highly important. Additionally, we perform some “stress tests” to check the robustness of the interaction of the different underlying methods. We begin with a repeating pattern of hourly rainfall sums of 0 mm and 40 mm and zero evaporation:

>>> inputs.rainfall.series = 12 * [0.0, 40.0]
>>> inputs.evaporation.series = 0.0
>>> inputs.capillaryrise.series = 0.0

The following “no evaporation” examples focus on a sufficiently deep soil column of 1 m depth:

>>> soildepth(1000.0)

three bins

For a series of 40 mm rainfall impulses, separated by hiatus periods of only one hour, one should expect infiltration rates to decrease from event to event. However, the eleventh event (21:00 - 22:00) deviates from this expectation. Inspecting the evolution of the wetting fronts’ relative moisture and depth reveals the cause. For the first event (01:00 - 02:00), a single saturated wetting front (handled by the second bin) is sufficient. For the next few subevents, an additional bin is required for modelling new wetting fronts until they reach the old redistributing wetting front of the second bin. From the seventh event (13:00 - 14:00), the new wetting front cannot reach the depth of the old wetting front during the event anymore, but still in the subsequent hiatus period. After the tenth event (19:00 - 20:00), the old front has sufficiently grown to be out of reach for the new front, even during the hiatus period. Then, ga_garto would need to activate a new bin to model the eleventh event with two old and one new wetting front, which is impossible due to the limited number of available bins:

ToDo: Should we think about increasing the number of bins automatically?

>>> test("ga_garto_24h_1000mm_3bins")
Click to see the table
Click to see the graph

There is no indication of an error in the water balance:

>>> round_(model.check_waterbalance(conditions))
0.0

four bins

After setting the number of available bins to four, ga_garto can create all wetting fronts required during the simulation period and predicts the infiltration rate decrease expected for all events:

>>> nmbbins(4)
>>> test("ga_garto_24h_1000mm_4bins")
Click to see the table
Click to see the graph

There is no indication of an error in the water balance:

>>> round_(model.check_waterbalance(conditions))
0.0

24 hours, evaporation

This subsection combines complex precipitation events with an evaporation rate of 1 mm/h:

>>> inputs.evaporation.series = 1.0

multiple events

In contrast to the deep soil example, rainfall starts in the second simulation step. Hence, the initial moisture (handled by the first bin) decreases during the first simulation step:

>>> test("ga_garto_24h_1000mm_evap")
Click to see the table
Click to see the graph

There is no indication of an error in the water balance:

>>> round_(model.check_waterbalance(conditions))
0.0

continuous event

Now, we define an extreme rainfall event with multiple increases and decreases in intensity:

>>> inputs.rainfall.series = [
...     0.0, 10.0, 20.0, 10.0, 20.0, 30.0, 20.0, 30.0, 40.0, 30.0, 40.0, 50.0,
...     50.0, 40.0, 30.0, 40.0, 30.0, 20.0, 30.0, 20.0, 10.0, 20.0, 10.0, 0.0]

Because of the vast rainfall, the soil saturates entirely in the 20th hour, and percolation becomes significant afterwards. In the 21st hour, it is restricted by saturated conductivity and in the 22nd hour, by the current rainfall rate. Note there is a difference in the behaviour especially relevant during the 21st hour, which we discuss in the documentation on method Redistribute_Front_V1:

>>> test("ga_garto_24h_1000mm_evap_continuous")
Click to see the table
Click to see the graph

There is no indication of an error in the water balance:

>>> round_(model.check_waterbalance(conditions))
0.0

long interruption

We define three strong rainfall impulses with hiatus periods of ten and eleven hours to give better insight into the longer-term behaviour of wetting front redistribution:

>>> inputs.rainfall.series = 0.0
>>> inputs.rainfall.series[[0, 11, 23]] = 40.0, 15.0, 15.0

The second, smaller rainfall impulse cause the creation of a second wetting front. This front does not become fully saturated and requires about three hours for to reach the first front’s depth:

>>> test("ga_garto_24h_1000mm_evap_interruptions")
Click to see the table
Click to see the graph

There is no indication of an error in the water balance:

>>> round_(model.check_waterbalance(conditions))
0.0

sealed surface

One can mark individual soil compartments as “sealed”:

>>> sealed(True)

Then, the soil’s surface prevents infiltration and all rainfall evaporates, or becomes surface runoff:

>>> test("ga_garto_24h_sealed")
Click to see the table
Click to see the graph

There is no indication of an error in the water balance:

>>> round_(model.check_waterbalance(conditions))
0.0

three compartments

We prepare three soil compartments within the same model instance, of which the last one is sealed:

>>> nmbsoils(3)
>>> sealed(False, False, True)

We must define each compartment’s area but only the depth of unsealed ones:

>>> soilarea(1.0, 2.0, 3.0)
>>> soildepth(400.0, 100.0, nan)

We take the soil properties for sand (first compartment) and clay (second compartment) from table 1 of Lai et al. (2015):

>>> residualmoisture(0.02, 0.09, nan)
>>> saturationmoisture(0.417, 0.385, nan)
>>> saturatedconductivity(235.6, 0.6, nan)
>>> poresizedistribution(0.694, 0.165, nan)
>>> airentrypotential(72.6, 373.0, nan)

Due to changing the arrays’ shapes, we need to update the memorised initial conditions required for checking the water balance:

>>> test.reset_inits()
>>> conditions = model.conditions

Sand allows for much more infiltration than clay, of course. So, the wetting front created during the first rainfall impulse reaches the sand soil’s bottom (at a depth of 400 mm) during the 6th hour. The tiny amount of percolation (0.0002 mm) at this time is due to numerical inaccuracy (or the fact that Merge_SoilDepthOvershootings_V1 considers water overshooting the soil’s depth as groundwater recharge). In contrast, in the clay soil, the deepest wetting front propagates to a depth of less than 70 mm during the entire 24 h period only:

>>> test("ga_garto_5h_3soils")
Click to see the table
Click to see the graph

There is no indication of an error in the water balance:

>>> round_(model.check_waterbalance(conditions))
0.0
class hydpy.models.ga_garto.Model[source]

Bases: BaseModel

HydPy-GA-GARTO (main model for calculating the Green-Ampt / Talbot-Ogden infiltration with redistribution).

The following “run methods” are called in the given sequence during each simulation step:
The following “additional methods” might be called by one or more of the other methods or are meant to be directly called by the user:
DOCNAME: DocName = ('GA-GARTO', 'main model for calculating the Green-Ampt / Talbot-Ogden infiltration with redistribution')
REUSABLE_METHODS: ClassVar[tuple[type[ReusableMethod], ...]] = ()
class hydpy.models.ga_garto.AideSequences(master: Sequences, cls_fastaccess: type[TypeFastAccess_co] | None = None, cymodel: CyModelProtocol | None = None)

Bases: AideSequences

Aide sequences of model ga_garto.

The following classes are selected:
class hydpy.models.ga_garto.ControlParameters(master: Parameters, cls_fastaccess: type[FastAccessParameter] | None = None, cymodel: CyModelProtocol | None = None)

Bases: SubParameters

Control parameters of model ga_garto.

The following classes are selected:
class hydpy.models.ga_garto.DerivedParameters(master: Parameters, cls_fastaccess: type[FastAccessParameter] | None = None, cymodel: CyModelProtocol | None = None)

Bases: SubParameters

Derived parameters of model ga_garto.

The following classes are selected:
class hydpy.models.ga_garto.FluxSequences(master: Sequences, cls_fastaccess: type[TypeFastAccess_co] | None = None, cymodel: CyModelProtocol | None = None)

Bases: FluxSequences

Flux sequences of model ga_garto.

The following classes are selected:
class hydpy.models.ga_garto.InputSequences(master: Sequences, cls_fastaccess: type[TypeFastAccess_co] | None = None, cymodel: CyModelProtocol | None = None)

Bases: InputSequences

Input sequences of model ga_garto.

The following classes are selected:
class hydpy.models.ga_garto.LogSequences(master: Sequences, cls_fastaccess: type[TypeFastAccess_co] | None = None, cymodel: CyModelProtocol | None = None)

Bases: LogSequences

Log sequences of model ga_garto.

The following classes are selected:
class hydpy.models.ga_garto.StateSequences(master: Sequences, cls_fastaccess: type[TypeFastAccess_co] | None = None, cymodel: CyModelProtocol | None = None)

Bases: StateSequences

State sequences of model ga_garto.

The following classes are selected:

  • Moisture() The relative soil moisture of each bin (within the wetting front) [-].

  • FrontDepth() The depth of the wetting front in each bin [-].