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.

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:

\theta_{WP}= \theta_{1500t} + (0.14 \theta_{1500t} - 0.002) with:

\theta_{1500t} = -0.024 S + 0.487 C + 0.006 OM + 0.005(S \times OM) - 0.013 (C \times OM) + 0.068 (S \times C) + 0.031

where:

  • S: represents the sand content of the soil (mass percentage),
  • C: represents the clay content of the soil (mass percentage),
  • OM: represents the organic matter content of the soil (mass percentage).
Similarly, the calculation of the water content at field capacity \theta_{FC} can be done as follows:
\theta_{FC} = \theta_{33t} + (1.283 \theta_{33t}^{2} - 0.374 \theta_{33t}-0.15) with:
\theta_{33t} = -0.251 S + 0.195 C + 0.011 OM + 0.006 (S \times OM) - 0.027 (C \times OM) + 0.452 (S \times C) + 0.299

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:

0M=1.724×OC

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:

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:

Looking at both pandas.DataFrame shows the following points:

  • the time resolution between both datasets is not the same,
  • for some reasons, potential evapotranspiration cannot be calculated at some dates. It corresponds to lines where the quality indicator ET_QC is higher than 1.

Both issues must be handled before implementing the iterative process: we want to work on a similar timeline with potential evapotranspiration and precipitation, and we want to avoid missing values.

Resampling the time resolution of an ee.ImageCollection

To address these issues (homogeneous time index and missing values), we make a sum resampling of both datasets by month. When PET cannot be calculated, the monthly averaged value is considered. The key steps and functions used to resample are described below:

  • A new date index is defined as a sequence using the ee.List.sequence method.
  • A function representing the resampling operation is defined. This function consists of grouping images of the desired time interval and calculating the sum. The sum is calculated by taking the mean between available images and multiplying it by the duration of the interval.
  • The user-supplied function is then mapped over the new time index using .map().

Finally, the resampling procedure reads as follows:

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:

time pr pet apwl st rech
datetime
2015-01-01 1420070400000 83.127947 23.637500 0.000000 2.914923e+01 59.490448
2015-02-01 1422748800000 61.235233 39.200000 0.000000 2.914923e+01 22.035233
2015-03-01 1425168000000 60.416391 76.434375 527.000626 3.043201e-24 0.000000
2015-04-01 1427846400000 48.404174 137.000000 447.422202 8.249401e-17 0.000000
2015-05-01 1430438400000 70.076297 187.550000 237.667646 1.357830e-04 0.000000
2015-06-01 1433116800000 89.217610 220.500000 131.282390 3.225894e-01 0.000000
2015-07-01 1435708800000 36.950253 228.915625 191.965372 4.022851e-02 0.000000
2015-08-01 1438387200000 56.904872 186.871875 321.932375 6.427792e-07 0.000000
2015-09-01 1441065600000 131.345901 108.656250 7.302447 2.268965e+01 0.000000
2015-10-01 1443657600000 84.328574 48.179167 0.000000 2.914923e+01 29.689833
2015-11-01 1446336000000 35.439609 39.500000 4.060391 2.535895e+01 0.000000
2015-12-01 1448928000000 16.603715 25.058333 12.515009 1.650707e+01 0.00000

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

Groundwater recharge comparison between multiple places

We now may want to get local information about groundwater recharge and/or map this variable on an area of interest.

Let’s define a function to get the local information based on the ee.ImageCollection we’ve just built:

We now use this function on a second point of interest located near the city of Montpellier (France). This city is located in the south of France, and precipitation and groundwater recharge are expected to be much lower than in the previous case.

pr pet apwl st rech
datetime
2015-12-31 602.685190 1479.557292 3469.362932 94.357806 22.609314
2016-12-31 668.897704 1449.505208 3515.844439 114.813998 86.200715
2017-12-31 371.787006 1538.359375 7834.083159 45.534037 16.356442
2018-12-31 826.732359 1458.902083 3520.025304 155.748515 147.778644
2019-12-31 550.755352 1454.181250 5354.379859 111.292106 121.659726

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.