This article is an Earth Engine community tutorial.

Introduction

Groundwater recharge represents the amount of water coming from precipitation reaching the groundwater table. Its determination helps to better understand the available/renewable groundwater in watersheds and the shape of groundwater flow systems.

One of the simplest methods to estimate groundwater recharge is the Thornthwaite-Mather procedure (Steenhuis and Van Der Molen, 1986). This procedure was published by Thornthwaite and Mather (1955, 1957). The idea of this procedure is to calculate the water balance in the root zone of the soil where water can be (1) evaporated into the atmosphere under the effect of heat, (2) transpired by vegetation, (3) stored by the soil, and eventually (4) infiltrated when stored water exceeds the field capacity.

This procedures relies on several parameters and variables described as follows:

• information about soil texture (e.g. sand and clay content) to describe the hydraulic properties of the soil and its capacity to store/infiltrate,
• meteorological records: precipitation and potential evapotranspiration.

Of course groundwater recharge can be influenced by many other factors such as the slope of the terrain, the snow cover, the variability of the crop/land cover and the irrigation. In the following these aspects are not taken into account.

In the first part of the tutorial, the Earth Engine python API will be initialized, some useful libraries will be imported, and the location/period of interest will be defined.

In the second part, OpenLandMap datasets related to soil properties will be explored. The wilting point and field capacity of the soil will be calculated by applying some mathematical expressions to multiple images.

In the third part, evapotranspiration and precipitation datasets will be imported. A function will be defined to resample the time resolution of an ee.ImageCollection and to homogenize time index of both datasets. Both datasets will then be combined into one.

In the fourth and final part, the Thornthwaite-Mather(TM) procedure will be implemented by iterating over the meteorological ee.ImageCollection. Finally, a comparison between groundwater recharge in two places will be described and the resulting mean annual groundwater recharge will be displayed over France.

Run me first

Earth Engine API

First of all, run the following cell to initialize the API. The output will contain instructions on how to grant this notebook access to Earth Engine using your account.

Other libraries

Import other libraries/modules used in this notebook.

• pandas: data analysis (including the DataFrame data structure)
• matplotlib: data visualization library
• numpy: array-processing package
• folium: interactive web map
• pprint: a pretty printer
• branca.colormap: utility module for dealing with colormaps.

Some input parameters

We additionally define some parameters used to evaluate the results of our implementation. Particularly:

• the period of interest to get meteorological records,
• a location of interest based on longitude and latitude coordinates. In the following, the point of interest is located in a productive agricultural region which is about 30 kilometers outside of the city of Lyon (France). This point is used to evaluate and illustrate the progress of the described procedure.

From soil texture to hydraulic properties

Definitions

Two hydraulic properties of soil are commonly used in the TM procedure:

• the wilting point represents the point below what water cannot be extracted by plant roots,
• the field capacity represents the point after which water cannot be stored by soil any more. After that point, gravitational forces become too high and water starts to infiltrate the lower levels.

Some equations given by Saxton & Rawls (2006) are used to link both parameters to the texture of the soil. The calculation of water content at wilting point \theta_{WP} can be done as follows:

We now define a function to get the ee.Image associated to the parameter we are interested in (e.g. sand, clay, organic carbon content, etc.).

We apply this function to import soil properties:

To illustrate the result, we define a new method for handing Earth Engine tiles and using it to display the clay content of the soil at a given reference depth, to a Leaflet map.

Now, a function is defined to get soil properties at a given location. The following function returns a dictionary indicating the value of the parameter of interest for each standard depth (in centimeter). This function uses the ee.Image.sample method to evaluate the ee.Image properties on the region of interest. The result is then transferred client-side using the ee.Image.getInfo method.

In the example below, we are asking for the sand content.

Sand content profile at the location of interest: {'b0': 0.36, 'b10': 0.35, 'b100': 0.37, 'b200': 0.39, 'b30': 0.34, 'b60': 0.35}

We now apply the function to plot the profile of the soil regarding sand and clay and organic carbon content at the location of interest:

Expression to calculate hydraulic properties

Now that soil properties are described, the water content at the field capacity and at the wilting point can be calculated according to the equation defined at the beginning of this section. Please note that in the equation of Saxton & Rawls (2006), the wilting point and field capacity are calculated using the Organic Matter content (OM) and not the Organic Carbon content (OC). In the following, we convert OC into OM using the corrective factor known as the Van Bemmelen factor:

Several operators are available to perform basic mathematical operations on image bands: add(), subtract(), multiply() and divide(). Here, we multiply the organic content by the Van Bemmelen factor. It is done using the ee.Image.multiply method on the organic carbon content ee.Image.

When the mathematical operation to apply to the ee.Image becomes too complex, the ee.Image.expression is a good alternative. We use it in the following code block since the calculation of wilting point and field capacity relies on multiple parameters and images. This method takes two arguments:

• a string formalizing the arithmetic expression we want to evaluate,
• a dict associating images to each parameter of the arithmetic expression.

The mathematical expression is applied as follows to determine wilting point and field capacity:

Let’s see the result around our location of interest:

Wilting point profile: {'b0': 0.152, 'b10': 0.158, 'b100': 0.175, 'b200': 0.169, 'b30': 0.175, 'b60': 0.181} 
Field capacity profile: {'b0': 0.271, 'b10': 0.278, 'b100': 0.289, 'b200': 0.28, 'b30': 0.294, 'b60': 0.297}

The result is displayed using barplots as follows:

Getting meteorological datasets

Datasets exploration

The meteorological data used in our implementation of the TM procedure relies on the following datasets:

• Climate Hazards Group InfraRed Precipitation with Station data (CHIRPS) gives precipitations on a daily basis (resolution of 5 km),
• MODIS Terra Net gives evapotranspiration on a 8-days basis (resolution of 500 m).

Both datasets are imported as follows, specifying the bands of interest using .select() and the period of interest using .filterDate().

Now we can have a closer look around our location of interest. To evaluate the properties of an ee.ImageCollection, the ee.ImageCollection.getRegion method is used and combined with ee.ImageCollection.getInfo method for a client-side visualization.

[['id', 'longitude', 'latitude', 'time', 'precipitation'], 
['20150101', 5.142855001584261, 45.77365530231022, 1420070400000, 0], 
['20150102', 5.142855001584261, 45.77365530231022, 1420156800000, 0], 
['20150103', 5.142855001584261, 45.77365530231022, 1420243200000, 0], 
['20150104', 5.142855001584261, 45.77365530231022, 1420329600000, 0]]

We now establish a procedure to get meteorological data around a given location in form of a pandas.DataFrame:

We apply the function and see the head of the resulting pandas.DataFrame:

We do the same for potential evaporation:

The precipitation dataset is now resampled by month as follows:

• the collection to resample is defined as pr,
• we want a collection on a monthly basis then freq = 1 and unit = "month",
• there is no correction factor to apply according to the dataset description then scale_factor = 1,
• "pr" is the name of the output band.

[['id', 'longitude', 'latitude', 'time', 'pr'], 
['0', 5.142855001584261, 45.77365530231022, 1420070400000, 83.12794733047485], 
['1', 5.142855001584261, 45.77365530231022, 1422748800000, 61.235233306884766], 
['2', 5.142855001584261, 45.77365530231022, 1425168000000, 60.416391015052795], 
['3', 5.142855001584261, 45.77365530231022, 1427846400000, 48.404173851013184]]

For evapotranspiration, we have to be careful with the unit. The dataset gives us an 8-day sum and a scale factor of 10 is applied. Then, to get a homogeneous unit, we need to rescale by dividing by 8 and 10:

[['id', 'longitude', 'latitude', 'time', 'pet'], 
['0', 5.142855001584261, 45.77365530231022, 1420070400000, 23.637500000000003], 
['1', 5.142855001584261, 45.77365530231022, 1422748800000, 39.2], 
['2', 5.142855001584261, 45.77365530231022, 1425168000000, 76.434375], 
['3', 5.142855001584261, 45.77365530231022, 1427846400000, 137]]

We now combine both ee.ImageCollection objects (pet_m and pr_m) using the ee.ImageCollection.combine method. Note that corresponding images in both ee.ImageCollection objects need to have the same time index before combining.

[['id', 'longitude', 'latitude', 'time', 'pr', 'pet'], 
['0', 5.142855001584261, 45.77365530231022, 1420070400000, 83.12794733047485, 23.637500000000003], 
['1', 5.142855001584261, 45.77365530231022, 1422748800000, 61.235233306884766, 39.2], 
['10', 5.142855001584261, 45.77365530231022, 1446336000000, 35.43960928916931, 39.5], 
['11', 5.142855001584261, 45.77365530231022, 1448928000000, 16.603714764118195, 25.058333333333337]]

We evaluate the result on our location of interest:

Implementation of the TM procedure

Description

Some additional definitions are needed to formalize the Thornthwaite-Mather procedure. The following definitions are given in accordance with Allen et al. (1998) (the document can be downloaded here):

TAW = 1000 \times (\theta_{FC} − \theta_{WP})\times Z{r} where:

• TAW: the total available soil water in the root zone (mm)
• \theta_{FC}: the water content at the field capacity (m^{3} m^{-3}),
• \theta_{WP}: the water content at wilting point (m^{3} m^{-3}),
• Z_{r}: the rooting depth (mm),

Typical values of \theta_{FC} and \theta_{WP} for different soil types are given in Table 19 of Allen et al. (1998).

The readily available water (RAW) is given by RAW = p \times TAW, where p is the average fraction of TAW that can be depleted from the root zone before moisture stress (ranging between 0 to 1). This quantity is also noted STFC which is the available water stored at field capacity in the root zone.

Ranges of maximum effective rooting depth Z_{r}, and soil water depletion fraction for no stress p, for common crops are given in the Table 22 of Allen et al. (1998). In addition, a global effective plant rooting depth dataset is provided by Yang et al. (2016) with a resolution of 0.5° by 0.5° (see the paper here and the dataset here).

According to this global dataset, the effective rooting depth around our region of interest (France) can reasonably assumed to Z_{r} = 0.5. Additionally, the parameter p is also assumed constant and equal to and p = 0.5 which is in line with common values described in Table 22 of Allen et al. (1998).

In the following, we also consider an averaged value between reference depths of the water content at wilting point and field capacity:

The Thornthwaite-Mather procedure used to estimate groundwater recharge is explicitly described by Steenhuis and Van der Molen (1985). This procedure uses monthly sums of potential evaporation, cumulative precipitation, and the moisture status of the soil which is calculated iteratively. The moisture status of the soils depends on the accumulated potential water loss (APWL). This parameter is calculated depending on whether the potential evaporation is greater than or less than the cumulative precipitation. The procedure reads as follow:

Case 1: potential evapotranspiration is higher than precipitation.

In that case, PET>P and APWL_{m} is incremented as follows: APWL_{m} = APWL_{m−1} + (PET_{m} − P{m})  where:

• APWL_{m} represents the accumulated potential water loss for the month m,
• PET_{m} the cumulative potential evapotranspiration at month m,
• P_{m} the cumulative precipitation at month m,

and the relationship between APWL and the amount of water stored in the root zone for the month m is expressed as: ST_{m} = STFC \times exp(−APWL_{m}/STFC) where ST_{m} is the available water stored in the root zone for the month m.

Case 2: potential evapotranspiration is lower than precipitation.

In that case, PET<P  and ST_{m} is incremented as follows: ST_{m} = ST_{m-1} + (P_{m} - PET_{m}) .

Case 2.1: the storage ST_{m} is higher than the water stored at the field capacity.

If ST_{m} > STFC the recharge is calculated as: R_{m}= ST_{m} + P_{m} - PET_{m}

In addition, the water stored at the end of the month m becomes equal to STFC and APWL_{m} is set equal to zero.

Case 2.2: the storage ST_{m} is less than or equal to the water stored at the field capacity.

If ST_{m} <= STFC, APWL_{m} is implemented as follows: APWL_{m} = STFC \times ln(ST_{m}/STFC), and no percolation occurs.

Initialization

The initial time of the calculation is defined according to the first date of the meteorological collection:

Then, we initialize the calculation with an ee.Image where all bands associated to the hydric state of the soil are set equal to ee.Image(0), except for the initial storage which is considered to be equal to the water content at field capacity, meaning that ST_{0} = STFC.

We combine all these bands into one ee.Image adding new bands to the first using the ee.Image.addBands method:

We also initialize a list in which new images will be added after each iteration. We create this server-side list using the ee.List method.

Iteration over an ee.ImageCollection

The procedure is implemented by means of the ee.ImageCollection.iterate method, which applies a user-supplied function to each element of a collection. For each time step, groundwater recharge is calculated using the recharge_calculator considering the previous hydric state of the soil and current meteorological conditions.

Of course, considering the TM description, several cases must be distinguished to calculate groundwater recharge. The distinction is made by the definition of binary layers with different logical operations. It allows specific calculations to be applied in areas where a given condition is true using the ee.Image.where method.

The function we apply to each element of the meteorological dataset to calculate groundwater recharge is defined as follows.

The TM procedure can now be applied to the meteorological ee.ImageCollection:

Let’s have a look at the result around the location of interest:

The result can be displayed in the form of a barplot as follows:

The result shows the distribution of precipitation, potential evapotranspiration, and groundwater recharge along the year. It shows that in our area of interest, groundwater recharge generally occurs from October to March. Even though a significant amount of precipitation occurs from April to September, evapotranspiration largely dominates because of high temperatures and sun exposure during these months. The result is that no percolation into aquifers occurs during this period.

Now the annual average recharge over the period of interest can be calculated. To do that, we resample the DataFrame we’ve just created:

The mean annual recharge at our point of interest is 154 mm/an

The result shows that the annual recharge in Lyon is almost twice as high as in the area of Montpellier. The result also shows a great variability of the annual recharge ranging from 98 mm/y to 258 mm/y in Lyon and from 16 mm/y to 147 mm/y in Montpellier.

Groundwater recharge map of France

To get a map of groundwater recharge around our region of interest, let’s create a mean composite ee.Image based on our resulting ee.ImageCollection.

And finally the map can be drawn.

The resulting map represents the averaged annual precipitation and groundwater recharge rate in France over the period 2015-2019. Groundwater recharge is represented between 0 (red) to 500 mm/year (purple). Precipitation is represented between 500 (white) to 1500 mm/year (blue).

Please note the following points:

• The precipitation dataset used in this tutorial is not available over +/- 50 deg latitudes. That is why (1) the northern part of France is not covered by the precipitation/recharge, and (2) some anomalies/discontinuities can be observed near the boundary of the dataset.
• The procedure described in this tutorial does not consider the terrain slope which has a strong influence on runoff, especially in areas with significant relief. In these areas, groundwater recharge is overestimated because runoff is not considered.
• Snow cover has a more complex contribution to groundwater recharge than precipitation. This process is not taken into account in this tutorial. Hence, groundwater recharge is overestimated in areas where snow cover period is significant.
• The recharge calculation can be highly influenced by missing precipitation or potential evapotranspiration values. This issue is partially addressed by the dataset resampling procedure, but missing value can still exist and this point must be checked before considering that the calculated recharge is representative and reliable.

Bibliography

Allen RG, Pereira LS, Raes D, Smith M (1998). Crop evapotranspiration: guidelines for computing crop water requirements. Irrigation and Drainage Paper 56, FAO, Rome.

Saxton, K. E., & Rawls, W. J. (2006). Soil water characteristic estimates by texture and organic matter for hydrologic solutions. Soil science society of America Journal, 70(5), 1569-1578.

Steenhuis, T. S., & Van der Molen, W. H. (1986). The Thornthwaite-Mather procedure as a simple engineering method to predict recharge. Journal of Hydrology, 84(3-4), 221-229.

Thornthwaite, C. W., & Mather, J. R. (1957). Instructions and tables for computing potential evapotranspiration and the water balance. Publ. Climatol., 10(3).

Yang, Y., Donohue, R. J., & McVicar, T. R. (2016). Global estimation of effective plant rooting depth: Implications for hydrological modeling. Water Resources Research, 52(10), 8260-8276.

Acknowledgements

Thanks to Susanne Benz for precious advice. Thanks to Justin Braaten and Gino Miceli for their help in the final production of this Earth Engine community tutorial.

Play with the script…

This tutorial is available as a notebook and can be run on Goolge Colab.