A simple method for water balance estimation based on the empirical method and remotely sensed evapotranspiration estimates

Estimation of the water balance is a complicated task since some of its constituents, especially changes in groundwater storage, are difficult to measure directly and are often estimated indirectly through various hydrological models or using empirical methods. Empirical methods have been developed for the computation of various water budget components, such as evapotranspiration (ET), groundwater recharge (GR) and surface runoff. A well-known and useful formula for surface runoff computation is the Soil Conservation Service Curve Number (SCS-CN) documented in Section 4 of the National Engineering Handbook (NEH-4) of the US Department of Agriculture in 1956 (Mishra & Singh 2003). A semi-empirical method for the estimation of potential monthly ET is the Thornthwaite and Mather's (T&M's) formula which appeared in 1955 and was revised in 1957 (Black 2008). Despite the severe criticism the method received, especially the 1955 edition which comprised some flaws, the specific approach based on the revisions of 1957 has been widely used with quite acceptable results (Black 2008). Nevertheless, T&M's formula does not account for vegetative effects. GR is the one component of the water budget which is almost impossible to be measured directly due to the nature of the property itself but also due to the heterogeneity and complexity of the subsurface (Kinzelbach et al. 2002). Although numerical modeling offers a robust tool for indirect estimation of GR, its applicability is often constrained due to data scarcity, as even where there are operational monitoring networks for hydrological parameters, the maintenance costs result in frequent failures to provide continuous hydrologic observations that would support the development of reliable hydrological models.

The need to obtain reliable GR estimates has been met by the development of several empirical methods which are usually applicable to the specific geographical regions where they are developed. Such methods usually estimate GR through a regression equation with precipitation and they offer a reasonable alternative for GR estimation especially in the United States and in India, such as the Maxey–Eakin method (Maxey & Eakin 1949), the Sehgal method (Sehgal 1973), the Hearne and Dewey (Hearne & Dewey 1988), the Anderson's formula (Anderson et al. 1992), the Waltemeyer approach (Waltemeyer 2001) and the methodology developed by Kambhammettu et al. (2011). A common feature of those methods is that they are developed for alluvial aquifers which are characterized by mountain system recharge. However, they implicitly assume that increased precipitation leads to increased GR, neglecting the impact of increased temperature which might transform excess rainfall to ET (Gemitzi et al. 2017). Those empirical methods have minimum data requirements, i.e. they usually require only precipitation data, but their verification needs either output from hydrological models or estimation of the rest of the water budget components, i.e. surface runoff and ET.

Advances in remote sensing have greatly contributed to modeling and prediction of basins with little or no monitoring, where continuous data records of stream discharge and groundwater levels are not available (Lakshmi 2013; Mohanty et al. 2013; Mohammed et al. 2018c, 2018a, 2018b). Previous research has shown that remotely sensed observations can be used to accurately compute the water balance by either incorporating remotely sensed precipitation data (Lu et al. 2018) or computing ET losses based on remotely sensed land surface temperature (Dalezios et al. 2018). Recently, Gemitzi et al. (2017) presented an empirical method for monthly GR computation, using a regression formula and monthly precipitation and remotely sensed monthly actual ET data. This latter work has been applied and verified in a catchment in NE Greece, demonstrating thus its applicability in a typical mountain front system recharge Mediterranean basin.

The present work is an extension of the previously developed GR estimation methodology of Gemitzi et al. (2017). Thus, we expanded our previous work to estimate not only GR but all water balance components at the monthly time step. This was achieved by scaling moderate resolution imaging spectroradiometer (MODIS) terrestrial ET product using monthly precipitation data as a containing parameter for ET. The application of the proposed method was demonstrated in seven basins in NE Greece, whereas performance assessment was conducted in one of the study catchments in comparison to the output from the widely used hydrological model, namely Soil and Water Assessment Tool (SWAT) (Arnold et al. 1998; Neitsch et al. 2011).

Computation of the water balance at the catchment scale requires measured or computed data of at least three of its constituents. The fourth can then be estimated solving Equation (1) for the unknown parameter. As precipitation is the most commonly measured parameter, measured precipitation data were used whereas our methodology focuses on the computation of the rest of the three components of Equation (1) at the monthly time step. Therefore, in this work, ET including soil water losses was determined using the MODIS remotely sensed ET product, whereas GR was computed using a previously defined empirical formula developed in the study area (Gemitzi et al. 2017). Surface runoff was then estimated as the remaining term of Equation (1). In the following paragraphs, the description of the process to estimate ET, surface runoff and changes in groundwater storage is presented.

ET is a crucial factor of the water balance as it routes precipitation back to the atmosphere, either directly as evaporation or indirectly through the transpiration process of plants. In tropical and mid-latitude countries, ET has been reported to account for more than 50% of precipitation (Leopoldo et al. 1995; European Academies Science Advisory Council (EASAC) 2010). Remote sensing has proved to be a convenient way to acquire estimates of ET globally. MODIS ET products offer global gridded data sets at 500 m spatial resolution (https://modis.gsfc.nasa.gov/data/dataprod/mod16.php). 8-day composites of ET and potential ET (PET) are available through the MOD16A2 (Terra satellite) and MYD16A2 (Aqua satellite).

The initial algorithm for computation of MODIS ET (Mu et al. 2007) used the Penman–Monteith ET formula (Monteith 1965) and MODIS land cover, albedo, Leaf Area Index (LAI), Fraction of Photosynthetically Active Radiation (FPAR) and Enhanced Vegetation together with daily meteorological reanalysis data from NASA's Global Modeling and Assimilation Office. The methodology was further improved to account for night-time ET components (Mu et al. 2011; Running et al. 2019). Since the algorithms for MODIS ET estimation incorporates soil water losses, MODIS ET was assumed to comprise both ET and L terms of Equation (1).

In the present work, the 8-day ET composites provided in the MOD16A2 data set were downloaded and processed from January 2009 to June 2019. Processing of ET data comprises filtering based on the Quality Control (QC) layer accompanying each ET layer and then aggregating of ET values throughout each month to compute the monthly ET pixel values. As the specific ET product is the result of a complex algorithm which makes use of many other remotely sensed parameters stated above, the QC layer contains information for the corresponding input LAI/FPAR (MYD15A2H) granule of the same 8-day composite period and not for ET. Therefore, to reduce the uncertainty of pixel-based ET values, a spatial average of all pixels in each catchment area examined was computed and used for further estimation of water balance at the basin scale, instead of proceeding with the computation of the water balance at the pixel level.

An encountered problem was that during summer months, ET_MODIS exceeded monthly precipitation. In this case, monthly ET is scaled to equal monthly precipitation and GR and SR are set to zero. In all other cases, i.e. when a surplus of precipitation is available, Equation (2) is used to compute GR. Monthly SR is then computed as the remaining quantity from Equation (1). The methodology is presented in the flow chart of Figure 1.

To test the accuracy of the produced results, a comparison against modeled results of the monthly water balance in one of the seven study basins was performed. Vosvozis river basin was selected as a testing site since for this specific catchment there is a previously developed, calibrated and verified hydrological model, i.e. the SWAT model (Pisinaras et al. 2014). The model was run with the same parametrization for the period from January 2013 to June 2019 and results for ET, GR and surface runoff were compared to the methodology presented herein. In the SWAT model, ET was calculated using the Penman–Monteith (Monteith 1965) method and SR with the SCS-CN method. Three performance metrics were computed, i.e. Index of Agreement (IoA) (Willmott 1982), Pearson's linear Correlation Coefficient (PCC) and Root-Mean-Squared Error (RMSE) given by:

• 1.
  IoA:

  

  (3)
• 2.
  PCC:

  

  (4)
• 3.
  RMSE:

  

  (5)

  where n is the number of months in the examined time series, is the observed value (SWAT model), is the predicted value from this method, and are the mean values of SWAT observations and predictions of this method, respectively. All those indices measure the performance of a model comparing the predicted and observed values, and they express the spread in (Groenendijk et al. 2014). In cases when significant errors in the model exist, the RMSE is difficult to use to objectively assess the agreement of modeled and observed values. Therefore, the dimensionless IoA is proposed as an alternative to express the agreement in a more direct way (Willmott 1982). IoA ranges from 0 (no agreement) to 1 (perfect agreement).

The presented methodology was applied in seven catchments in the Thrace region (NE Greece) (Figure 2). The examined basins are (from west to east) Kosinthos, Vosvozis, Filiouris, Irini, Erithropotamos, South Evros and North Evros river basins and they were selected based on the precipitation data availability within or close to each basin. Mean annual precipitation (2013–2018) in Thrace ranges from 522 mm (Alexandroupolis meteo-station) to 667 mm (Didimoticho meteo-station). The climate of the area is Mediterranean with dry hot summers and wet winters (Gemitzi et al. 2017). The topography of the area comprises a flat terrain in the southern and eastern parts and a mountainous area, i.e. Rhodope mountain range, in the northern parts and along the Greek–Bulgarian borders. Rhodope mountain range is a pristine area covered by broad-leaved and coniferous forest. Rhodope mountain range is formed by Paleozoic metamorphic rocks, consisting of gneisses, schists, amphibolites and marbles (Kilias et al. 1999). Fractured aquifers are formed within those hard rock formations. The plain areas are occupied by extended agricultural fields with the main urban centers of the region, i.e. Xanthi, Komotini, Alexandroupolis, Didimoticho and Orestias, also located in this area. Those flat regions are formed by alluvial Pliocene and Quaternary deposits, and porous aquifers are formed within those formations. Previous research (Gemitzi 2012; Gemitzi et al. 2017) has shown that the study area is a typical mountain system recharge region where infiltration of precipitation through mountain bedrock and stream bed infiltration of mountain system runoff (Ajami et al. 2012) are the dominating hydrological processes that recharge the alluvial aquifers in the plain terrain.

Remotely sensed MODIS 8-day ET (MOD16A2) were obtained for a 10-year period, i.e. 2009–2019 and were processed and aggregated at the monthly time step, whereas monthly precipitation data acquired from seven meteorological stations operating in the study region were also obtained (Figure 2). Not all meteo-stations operated for the same time period. Therefore, the period of available precipitation data for each one of the seven stations is shown in Table 1. All data for basin scale monthly MODIS ET values and monthly precipitation for each one of the seven monitoring sites are provided in the Supplementary Material accompanying the present work.

Since there is an uneven spatial distribution of precipitation data, spatial averages of monthly precipitation data were acquired for each study basin aggregating precipitation as monthly weighted averages of point data values using Thiessen polygons (Thiessen 1911). Spatial interpolation with Thiessen polygons is based on the assumption that each measuring station is surrounded by its area of influence where for any point within this area, precipitation is considered equal to the observed precipitation at the closest gauge (Schumann 1998). The main disadvantage of the method is that it does not account for topographic gradients; however, for relatively flat catchments, it is considered as a fairly accurate interpolation procedure (Olsson et al. 2013). In the case of the study catchments where the topography is not flat through the basins, a comparison against remotely sensed GPM remotely sensed monthly precipitation data of 0.1° spatial resolution (Huffman et al. 2019) was performed. Monthly GPM precipitation gridded data were averaged at the catchment scale and were compared to the Thiessen interpolated precipitation values for each examined catchment.

All computations were performed using the open R programing platform and its packages Raster ver. 2.9–23 (Hijmans 2017) for reading, filtering and aggregating at the basin scale MODIS ET data, deldir ver. 0.1–23 (Turner 2019) for creating spatial aggregates of precipitation data using Thiessen polygons and ggplot2 ver. 3.2.0 (Wickham et al. 2019) for preparing all graphics.

Comparison of the interpolated precipitation data with GPM data indicated that in all basins, there was a significant agreement of the two data sets as evidenced by the performance metrics examined (please refer to Supplementary Material). Therefore, IoA ranged from 0.997 (North Evros) to 0.975 (South Evros), PCC ranged from 0.993 (Vosvozis) to 0.954 (South Evros). RMSE ranged from 4.5 mm/month in North Evros to 11.2 mm/month in South Evros. Thus, the interpolation of in situ data does not seem to have introduced substantial errors in the computational process. Results of the methodology for the estimation of water balance components are presented in Figure 3. Since the operation period is different for each meteo-station, the computation period for each study basin is different. A parallel trend with different magnitudes is observed for both GR and SR. This is an expected outcome since Equation (2) partitions approximately 52% of effective precipitation into GR with the remaining part contributing to SR.

A common outcome observed in all basins is that ET is the largest component of the water budget, returning thus most of the precipitation back to the atmosphere. The lowest ET rates are observed in Vosvozis and Filiouris basins which might be associated with the lower precipitation rates in combination with the plain topography and alluvial deposits which favor GR. Results of the mean annual water balance in all seven basins for the period from 2013 to 2018 (common period of water balance computations for all study basins) are shown on the box plots of Figure 4 and can be found in tabular form in Supplementary Material.

According to these annual estimates of the water balance (Figure 4), ET in the study regions accounts for approximately 54% of precipitation (mean ET of all examined basins) during 2013–2018 (excluding 2019 as the complete precipitation record for this year was not available at the time of completion of this study). Mean annual ET in the study basins varied in the range of 39.3–62.5% of precipitation for the same period. GR was found to account for 24% of precipitation, ranging from 19.5% to 31.5% of precipitation, while SR has a mean annual precipitation of 22%, varying from 17.9% to 29.1% of precipitation. Those results are close to the ones estimated in National Program for the Development and Protection of Water Resources (Koutsoyiannis et al. 2008) and are published in the country report of the European Academies Science Advisory Council (EASAC) (2010). In that work, ET accounts for approximately 51% of precipitation and the sum of GR and SR accounts for remaining 49% of precipitation. Examining previously published research dealing with the estimation of water balance in areas all over Greece, the mean modeled values for runoff in 14 catchments in Thessaly (Central Greece) from 1960 to 2002 are reported to range from 54% to 15% of precipitation (Vasiliades et al. 2011). Mean values of water balance components for Aison River Basin (Central Greece) from 1976 to 2007 were reported in Karpouzos et al. (2011), indicating a mean actual ET of approximately 62%, and mean annual surface runoff and infiltration of 9% and 29%, respectively, which are comparable to the results found in the present work taking into account that all the above-reported areas are located south of the ones analyzed herein and at a different hydrogeological setting. In Gemitzi et al. (2017), GR for Vosvozis catchment was determined using environmental stable isotopes during summer 2013, and its values were reported to range from 10.7% to 31.9% of precipitation.

Performance of the methodology is evaluated through a comparison of the acquired results to SWAT model results, in the testing basin, i.e. Vosvozis basin, from January 2013 to June 2019. Figure 5 presents the comparison of ET, GR and SR computed with the present methodology and by SWAT model for the predefined 5-year period. Performance metrics are shown in Table 2.

Results shown in Figure 5 indicate a satisfactory performance of the proposed methodology which is also evidenced by the metrics presented in Table 2. In general, MODIS ET is higher than SWAT ET and ET scaled. ET scaled and SWAT ET demonstrate obvious similar trends, with SWAT ET demonstrating always higher values in dry periods compared to ET scaled (Figure 5(a)). This is because SWAT allows for ET from soil water even when precipitation is limited. The present methodology does not account for soil water ET since this is implicitly incorporated in the algorithm of MODIS ET computations (Mu et al. 2011). Scaling, however, MODIS ET with local precipitation data might result in the underestimation of soil water losses during summer months as it forces losses to equal precipitation during those months. A large deviation of MODIS ET from SWAT ET and ET scaled is observed during the summer months of 2016 and 2017, when there was very low precipitation in the study basin. In this case, MODIS ET seems to be closer to potential ET and not the actual. This may be attributed to the use of precipitation as a constraining parameter for actual ET computations in SWAT algorithm and in our methodology, which is not used in MODIS ET algorithm. Nevertheless, the accuracy of the acquired results (Table 2) indicates that scaling MODIS ET is a reasonable approximation for ET computation.

Concerning GR, there is an agreement between the SWAT model and our method (Figure 5(b)). SR is found to be higher in dry periods in the SWAT model, which corresponds to the base flow (Figure 5(c)). However, the trends of the two data sets also this time are similar. Performance metrics in Table 2 indicate a very high IoA meaning that the time series of water balance components evaluated with SWAT model and with the present methodology demonstrate common trends. Furthermore, the two other metrics, i.e. PCC and RMSE, indicate satisfactory outcomes. A considerably lower performance was detected for the raw (unscaled) ET MODIS data, due to the higher ET values, especially during summer, compared to local precipitation data. This finding can be explained examining the MODIS ET computation algorithm (Mu et al. 2011; Running et al. 2019), that makes use of NASA's Global Modeling and Assimilation Office reanalysis meteorological data with a coarse spatial resolution, which are afterwards interpolated to the spatial resolution of MOD16 products, i.e. 0.5 km. Even after being interpolated, such reanalysis data cannot describe the spatial variability of meteorological parameters at the local scale. Moreover, remotely sensed ET is estimated from data representing clear sky days only and there is no input from cloudy days and this is certainly a limitation of remotely sensed ET products. A recent validation study of MODIS ET (Running et al. 2019) indicated a mean absolute error ranging from 24.1% to 24.6%. Therefore, scaling remotely sensed ET with in situ meteorological data is recommended when local or regional research is the case. In our case, the scaling of MODIS ET was based only on precipitation information from stations operating in the plain part of the study area, as there is a lack of available meteorological information in the mountainous parts. Better coverage of the meteorological stations of the study area, especially of its mountain parts, would certainly improve the accuracy of the outcome.

The results of the methodology presented herein indicated that it can be used for computation of the water balance at the monthly or coarser time step at the basin scale. It is certain that it cannot depict the fine spatial and temporal resolution of the results of a hydrologic model, but it can accurately allocate monthly precipitation to various water components at the catchment scale with a reasonable accuracy. A limitation is that it was verified against a hydrologic model to only one basin. This was due to data scarcity in surrounding catchments for the development and calibration of a hydrologic model. A future goal therefore should be the application of the methodology to other Mediterranean areas where the comparison with hydrologic modeling results is feasible. An advantage of our method is that it requires only monthly precipitation information, either from earth-based meteorological stations or as remote sensing product. It can be applied even to ungauged basins, provided that there is an established relationship between precipitation and GR (usually from the previous application of a hydrological model) in a basin of analogous hydrological and climate setting.

In the present work, a simple methodology for the estimation of the individual monthly water balance components at the basin scale is proposed. The methodology is based on a previously developed empirical relationship between GR and effective precipitation. Remotely sensed MODIS ET data were processed and scaled before they were incorporated in the process, along with available monthly precipitation data. The proposed approach was applied in seven mountain system recharge basins in NE Greece for a 10-year period, i.e. 2009–2019, and it was verified in one basin by comparing its results with those acquired from the SWAT hydrological model. Three performance metrics were examined, which indicated that results are of reasonable accuracy. The individual monthly water balance components computed with this new approach at the basin scale resulted in mean ET losses accounting for approximately 54% of precipitation, mean GR of 24% and mean surface runoff of approximately 22% of precipitation in the study area. Although the methodology presented herein cannot depict the fine spatial and temporal resolution of hydrological models, it certainly offers a reliable and inexpensive alternative to hydrological modeling for monthly water balance computation at the catchment level, when little or no gauging stations exist.

The second author acknowledges the Technical Chamber of Greece for its support through the project ‘Climate Changes and Dependent Eco systems in Eastern Macedonia and Thrace’ (Project ID: 81472 – Democritus University of Thrace). We would like to thank the anonymous reviewers for their time and effort to review and improve our work.

The Supplementary Material for this paper is available online at https://dx.doi.org/10.2166/hydro.2020.182.