subsurface flow models

Journal of Hydrology xxx (2015) xxx–xxx

Contents lists available at ScienceDirect

Journal of Hydrology journal homepage: www.elsevier.com/locate/jhydrol

A simple iterative method for estimating evapotranspiration with integrated surface/subsurface flow models H.-T. Hwang a,⇑, Y.-J. Park a, S.K. Frey a, S.J. Berg a, E.A. Sudicky a,b a b

Aquanty, Inc., 564 Weber Street North, Unit 12, Waterloo, Ontario N2L 5C6, Canada Department of Earth and Environmental Sciences, University of Waterloo, 200 University Avenue West, Waterloo, Ontario N2L 3G1, Canada

a r t i c l e

i n f o

Article history: Received 3 April 2015 Received in revised form 14 July 2015 Accepted 3 October 2015 Available online xxxx This manuscript was handled by Konstantine P. Georgakakos, Editor-in-Chief, with the assistance of Jennifer Guohong Duan, Associate Editor Keywords: Evapotranspiration Integrated surface/subsurface model Stream flow Iterative method Precipitation HydroGeoSphere

s u m m a r y This work presents an iterative, water balance based approach to estimate actual evapotranspiration (ET) with integrated surface/subsurface flow models. Traditionally, groundwater level fluctuation methods have been widely accepted and used for estimating ET and net groundwater recharge; however, in watersheds where interactions between surface and subsurface flow regimes are highly dynamic, the traditional method may be overly simplistic. Here, an innovative methodology is derived and demonstrated for using the water balance equation in conjunction with a fully-integrated surface and subsurface hydrologic model (HydroGeoSphere) in order to estimate ET at watershed and sub-watershed scales. The method invokes a simple and robust iterative numerical solution. For the proof of concept demonstrations, the method is used to estimate ET for a simple synthetic watershed and then for a real, highlycharacterized 7000 km2 watershed in Southern Ontario, Canada (Grand River Watershed). The results for the Grand River Watershed show that with three to five iterations, the solution converges to a result where there is less than 1% relative error in stream flow calibration at 16 stream gauging stations. The spatially-averaged ET estimated using the iterative method shows a high level of agreement (R2 = 0.99) with that from a benchmark case simulated with an ET model embedded directly in HydroGeoSphere. The new approach presented here is applicable to any watershed that is suited for integrated surface water/groundwater flow modelling and where spatially-averaged ET estimates are useful for calibrating modelled stream discharge. Ó 2015 Elsevier B.V. All rights reserved.

1. Introduction The dynamics of water movement and distribution within the hydrologic cycle involves physical processes in the atmosphere, land surface and subsurface and their interactions (VanderKwaak and Loague, 2001; Panday and Huyakorn, 2004; Kollet and Maxwell, 2006; Kundzewicz et al., 2007; Maxwell et al., 2007; Aquanty Inc., 2015; Delfs et al., 2013; Niu et al., 2014). The processes controlling the movement of water from one reservoir to another include precipitation, evapotranspiration (ET), overland flow, infiltration, recharge/discharge, and groundwater flow. Within the water cycle, understanding the balance of available water between reservoirs is critical for analyzing the sustainability of the surface and subsurface water resources (Partington et al., 2011; Rassam et al., 2013; Li et al., 2015). The concept of a watershed provides a convenient logical unit for hydrologic analyses, as it serves as semi-closed system for ⇑ Corresponding author. Tel.: +1 519 279 1080x123; fax: +1 519 279 1081. E-mail address: [email protected] (H.-T. Hwang).

water. In an ideal watershed, precipitation is the only source of water; and thus, it is balanced by all of the sinks in the system, such as stream flow at the watershed outlet, ET, and anthropogenic water consumption for urban and agricultural purposes (Li et al., 2008; Bolger et al., 2011; Pérez et al., 2011; Frey et al., 2013; Condon and Maxwell, 2014; De Schepper et al., 2015). When average precipitation and stream discharge at the outlet are known for a natural watershed, ET can be accurately estimated if the assumption is made that groundwater flow is limited between the watershed and its surroundings. At the sub-watershed scale, the use of more sophisticated tools such as numerical modelling and isotope tracers are often required for water balance analysis, as these areas are typically considered as open groundwater flow systems (Yi et al., 2010; Bearup et al., 2014; Ala-aho et al., 2015). In analyzing the water balance for a hydrologic system, precipitation and stream flow are relatively easy to measure compared to groundwater flow, infiltration and exfiltration, and ET. Because changes in ET impact on the amount of soil water availability, ET is an important hydrologic component that controls water cycling in the eco-hydrosphere as well as the moisture distribution in the

http://dx.doi.org/10.1016/j.jhydrol.2015.10.003 0022-1694/Ó 2015 Elsevier B.V. All rights reserved.

Please cite this article in press as: Hwang, H.-T., et al. A simple iterative method for estimating evapotranspiration with integrated surface/subsurface flow models. J. Hydrol. (2015), http://dx.doi.org/10.1016/j.jhydrol.2015.10.003

2

H.-T. Hwang et al. / Journal of Hydrology xxx (2015) xxx–xxx

atmosphere (Sellers et al., 1996; Gowda et al., 2008; Hanson, 2009). Many variations of the water budget method have been used to estimate ET for watershed and river basin scale analyses (Quinn and Beven, 1993; Hornberger, 1998; Rodell and Famiglietti, 2002; Rodell et al., 2004; Buerge et al., 2009; Davis and Dukes, 2010; Van Stempvoort et al., 2011; Dastorani and Poormohammadi, 2012; Robertson et al., 2013; Xue et al., 2013; Hassan et al., 2014; Ala-aho et al., 2015). Integrated surface and subsurface hydrological models are becoming increasingly popular tools to characterize major hydrological processes, as well as detailed surface water and groundwater interactions within watersheds. The physical processes considered and the degree of comprehensiveness reflected in the implementation of those processes can vary widely among the different models; Ebel and Loague (2006) indicate that physics-based integrated surface/subsurface simulations can increase both the model complexity and the uncertainty associated with the simulation results on account of the additional parameters and calibration data that are required. For this reason, many studies have combined integrated numerical models with neural network and inverse algorithms to analyze hydrological processes (Maneta et al., 2008; Goderniaux et al., 2009). In this study, the main objectives are (1) to suggest a simple Newton iterative approach that estimates actual ET from measured stream flow using the HydroGeoSphere (HGS) integrated surface/subsurface model (Aquanty Inc., 2015; Hwang et al., 2014); and (2) to validate the accuracy and applicability of the iterative model with a realworld watershed.

2. Theory This study presents an iterative method to estimate ET using an integrated surface/subsurface model. The ET estimation is based on the water balance in a watershed system. Fig. 1 depicts the major processes controlling the water balance in an idealized watershed where two upstream and downstream sub-watersheds are assumed to interact with each other. The outflow of surface water (surfout) and groundwater (GWout) from the upstream subwatershed are the inflow for surface water (surfin) and groundwater (GWin) for the downstream sub-watershed. For the entire watershed and each of the sub-watersheds, precipitation (P) acts as a source of water, evapotranspiration (ET) is a sink, and the

surface and subsurface flow regimes interact through infiltration and exfiltration. Water balance at steady-state conditions in an integrated surface and subsurface flow system can be described by the following equation: in out out Q p þ Q in surf þ Q GW ¼ Q ET þ Q surf þ Q GW

ð1Þ

where the left side of the equation represents the sources of water for the system in terms of the average rate of effective precipitation (Qp) consisting of liquid precipitation and snowmelt, and surface in water and groundwater flowing into the domain ðQ in surf and Q GW ), and the right terms are the sinks and include ET (QET) and surface out water and groundwater flowing out of the system ðQ out surf and Q GW ). In Eq. (1), no storage change is considered as the system is assumed to be at equilibrium. Using an integrated surface/subsurface model with ET, the water balance can be presented in the model as

b in ðQ Þ þ Q b in ðQ Þ ¼ Q b ET ðQ Þ þ Q b out ðQ Þ þ Q b out ðQ Þ Qp þ Q p p p p p surf GW surf GW

ð2Þ

b represents the simulated Q using the numerical model. If where Q the stream flow rates are measured at specific locations in the system, an inverse problem can be defined in order to find the optimal parameterization for flow and ET that minimizes the difference between measured and simulated flow rates. It is noted that among the hydrologic processes simulated by physics-based integrated numerical models, ET typically requires the largest number of parameters as it is a function of both surface and subsurface moisture conditions as well as climate conditions (potential ET); and thus, this inverse problem often becomes an ill-posed problem. As this study is concerned with the average rate of ET, the use of net precipitation (Qp  QET) can significantly simplify the system. If the net precipitation (Qp  QET) is represented by aQp or a  (Qp  QET)/Qp, Eq. (2) can be re-arranged as:

b out aQ p þ Qb insurf ðaÞ þ Qb inGW ðaÞ ¼ Qb out surf ðaÞ þ Q GW ðaÞ

ð3Þ

where a is the ratio of net precipitation to total precipitation. In Eq. (3), surface water and groundwater flow components are assumed to be driven by a given net precipitation in the model. Generally, stream flow at a measurement location can be considered within the context of the water balance for the catchment area of that specific location. Taking into account the water balance at the stream flow measurement location (i.e. Q in surf ¼ 0 and

Fig. 1. Schematic of water balance in upstream and downstream watersheds.

Please cite this article in press as: Hwang, H.-T., et al. A simple iterative method for estimating evapotranspiration with integrated surface/subsurface flow models. J. Hydrol. (2015), http://dx.doi.org/10.1016/j.jhydrol.2015.10.003

3

H.-T. Hwang et al. / Journal of Hydrology xxx (2015) xxx–xxx

Q in GW ¼ 0), the relationship between the stream flow and the water balance can be rewritten as:

b out ðaÞ b out ðaÞ ¼ aQ  Q Q p surf GW

ð4Þ

b out ðaÞ is the simulated surface (or stream) flow rate at the where Q surf b out ðaÞ is the simulated groundwater measurement location, and Q GW discharge rate at the same location. The difference between the measured and simulated stream flow can therefore be defined as:

h i out b out b out f ðaÞ ¼ Q out surf  Q surf ¼ Q surf  aQ p  Q GW ðaÞ

ð5Þ

where f(a) represents the difference between simulated and measured surface flow rates. As the change in net precipitation ratio (Da) affects surface water and groundwater flow in the system, the derivative of Eq. (5) can be 0

f ðaÞ ¼ Q p þ

b out @Q GW @a

ð6Þ

A simple Newton iterative procedure can then be applied to estimate the average upstream net precipitation (a) as: 0

f ðak þ Dak Þ  f ðak Þ þ f ðak ÞDak  0

ð7aÞ

0

f ðak ÞDak ¼ f ðak Þ

ð7bÞ

akþ1 ¼ ak þ Dak

ð7cÞ

where k and k + 1 represent the current and next iteration levels, respectively. When the groundwater flow is significantly less responsive to net precipitation compared to the surface water flow, a further approximation can be applied for the derivative as

Qp 

b out @Q 0 GW and f ðaÞ  Q p @a

ð8Þ

The Newton iterative method for estimating a then can be simplified as

akþ1 ¼ ak þ

 k b out Q out surf  Q surf

ð9Þ

Qp

In cases where the stream flow is measured at multiple locations in a watershed (Q out i;surf , i ¼ 1 to nobs ), the average net precipitation, ai Q i;p , is calculated recursively for the contributing area represented by each measurement location such that:

akþ1 ¼ aki þ Daki i



Daki ¼

b out Q out i;surf  Q i;surf

k



ð10aÞ

  k  out b out j¼1 Q j;surf  Q j;surf

Pni

ð10bÞ

Q i;p

where Qi,p is the total precipitation over the catchment area for the stream flow measurement location, excluding the catchment area b out represent for any upstream measurement locations. Q out and Q i;surf

i;surf

the measured and simulated stream flows at location i, respectively. b out are the measured and simulated stream Similarly, Q out and Q j;surf

j;surf

flows at upstream location j, respectively, and ni is the number of further upstream measurements which affect Q out i;surf . Additionally, downstream flow measurements need to take into account the upstream measurements because areas affecting the downstream measurements overlap with those affecting to the upstream measurements. Therefore, the primary purpose of ni is to remove the effects of upstream flow. 3. Proof-of-concept demonstration Before demonstrating a real-world application, a simple conceptual watershed model is presented to validate the iterative

Fig. 2. Schematic of a conceptual watershed model.

method. The virtual watershed is 100 km long and 10 km wide and it is drained by a single stream located along the middle of the domain (Fig. 2). Land surface slopes parallel to the stream follow a quadratic function. The hydrologic scenario is based on conditions in Southern Ontario, Canada, with annual average precipitation of 900 mm/year and annual average potential evaporation of 700 mm/year assigned as flux boundary conditions to the surface of the model for the entire watershed. A control-volume finite element integrated surface/subsurface numerical model was constructed for the virtual watershed with the 2D overland domain represented by rectangular elements (200 m  200 m) and the subsurface domain discretized using hexahedral elements with variable thickness from 0.2 m (near surface) to 200 m (deep subsurface). The virtual watershed was divided into four subwatersheds (SW1 to SW4 from upstream to downstream) based on an assumption that four hypothetical stream flow measurement locations are uniformly spaced along the stream. Table 1 summarizes the flow and ET properties used for the reference simulation.

Table 1 Flow and ET properties for the reference simulation. Parameter

Value

Subsurface domain Kx = Ky = 10Kz (m/s)

5  106

Surface domain Roughness coefficient (m1/3 s) Coupling length (m)

0.095 102

Evapotranspiration Evaporation depth (m) Root depth (m) LAI table (–)

0.5 1.0 2.6

Transpiration fitting parameters (–) C1 C2 C3

0.12 0.2 2.0

Transpiration limiting saturation (–) Wilting point Field capacity Oxic limit Anoxic limit

0.3 0.65 0.9 1.0

Evaporation limiting saturation (–) Min Max

0.3 1.0

Canopy storage (–) Initial interception (–)

0.0 0.0

Please cite this article in press as: Hwang, H.-T., et al. A simple iterative method for estimating evapotranspiration with integrated surface/subsurface flow models. J. Hydrol. (2015), http://dx.doi.org/10.1016/j.jhydrol.2015.10.003

4

H.-T. Hwang et al. / Journal of Hydrology xxx (2015) xxx–xxx

Note that an ET and interception model suggested by Kristensen and Jensen (1975) is used in this study for computing the actual ET. Details of the ET model can be found in the HGS documentation (Aquanty Inc., 2015). For the validation of the iterative method suggested in this study, an integrated surface/subsurface flow simulation that includes ET processes was performed first to obtain simulated surface flow and ET values, and then the simulated surface flow rates at the four sub-watersheds were used as a calibration target for the surface/subsurface model with the iterative method, which estimates ET without the ET model. For the simulation with the iterative approach, all the parameter values other than the ET model parameters are the same as those of the reference case.

Fig. 3 shows the results of the reference simulation: spatial distributions of surface water depth in the overland flow domain, the rate of fluid exchange between surface and subsurface, the rate of actual ET, and the saturation and head in the subsurface domain. Regarding the estimation of average ET rates using the iterative method, if net precipitation is 83% of the precipitation for the entire watershed (a0i ¼ 0:83; i = 1–4), the stream flow estimated from the iterative method becomes the same as the simulated (reference model) stream flow at the outlet of the watershed, however, it is actually less than the simulated (reference model) stream flows at three upstream gauging stations by 37–3%. When the net precipitation applied to each sub-watershed is updated by Eq. (10), the difference between the estimated and simulated

Fig. 3. Results of integrated surface/subsurface simulation.

Please cite this article in press as: Hwang, H.-T., et al. A simple iterative method for estimating evapotranspiration with integrated surface/subsurface flow models. J. Hydrol. (2015), http://dx.doi.org/10.1016/j.jhydrol.2015.10.003

5

H.-T. Hwang et al. / Journal of Hydrology xxx (2015) xxx–xxx

to 24% (SW4) of the precipitation and the ratio between the stream flow at the outlet of each sub-watershed and the net precipitation applied to the sub-watershed ranges from 0.24 (SW1) to 1.53 (SW2). Fig. 4b shows that the estimated average rates of ET for each sub-watershed ranges from 52 mm/year (SW1) to 216 mm/year (SW4): the absolute estimation error is less than 2.5 mm/year and the relative error is about 3.8%, 1.5%, 0.9%, and 0.6% for each of the sub-watersheds, respectively. The error is relatively small and negligible compared to other uncertainties involved in the water budget analysis; such as ET parameters, precipitation, and heterogeneities of geological settings. Additionally, it should be noted that the estimated ET at SW4 is about four times higher than that at SW1 because the groundwater table is deeper and the soil saturation near ground surface is lower in the upstream subwatersheds (i.e. SW1 and SW2) than those in the downstream ones (i.e. SW3 and SW4). 4. Results and discussion 4.1. Grand River Watershed integrated hydrologic model

Fig. 4. Results comparison between simulated and estimated (a) surface flow rates and (b) ET at each subwatershed.

(reference model) stream flow rates at the uppermost gauging station is reduced to 9% and it becomes less than 0.1% at the rest of the downstream stations (Fig. 4a). A total of five iterations were required for the maximum difference to be less than 0.1% at all the stations. Table 2 summarizes the water budget for each of the subwatersheds as well as for the entire watershed based on the reference simulation. From the water budget for the entire watershed provided in Table 2, the average ET is calculated to be 17% of the precipitation while the outflow at the outlet accounts for 83% of precipitation. For the sub-watersheds, ET accounts for 6 (SW1) Table 2 Summary of water budget for the entire watershed. Parameter (m3/s)

Precipitation Evapotranspiration Surface flow

Sub-watershed

Total watershed

1

2

3

4

7.14 0.43 1.62

7.14 1.31 10.54

7.14 1.35 17.79

7.14 1.70 23.73

28.56 4.79 –

The iterative ET estimation approach proposed in this work is applied to a real world case study to demonstrate its applicability and robustness. The Grand River Watershed (GRW) in Southern Ontario, Canada is selected for the case study because the area has been well characterized for surface/subsurface hydrologic parameters including stream flow and hydrogeological unit information. Surface water and groundwater resources are heavily dependent upon within this 7000 km2 mixed use (agriculture, municipal, industrial) watershed, and within it, almost one million people are reliant upon groundwater for drinking water, which has provided a major impetus for detailed hydrologic characterization. For the Grand River Watershed, a numerical groundwater-only model was developed and calibrated by the Grand River Conservation Authority (AquaResource Inc., 2009). The model includes three major unconsolidated materials near the ground surface, and the underlying bedrock consists of 11 sedimentary formations, for which descriptions and hydraulic parameters are given in Table 3. The elevations of the bedrock surface range from 175 m in the south to 525 m in the north (Lake Erie Source Protection Region, 2008). Above the bedrock, the top seven layers of the model consist of surficial geology. For the integrated hydrologic simulations required for this demonstration, an HGS surface water/variablysaturated subsurface flow model was constructed based on the available 25 m resolution Digital Elevation Model (DEM) and land cover classification dataset (Grand River Conservation Authority, 1999), which were combined with the original subsurface model. The land cover dataset is based on remote imagery (Grand River

Table 3 Parameter values for the GRW subsurface domain (AquaResource Inc., 2009). Subsurface material type

Hydraulic conductivity Kx = Ky = 10Kz (m/s)

Till (top two element layers) Sand/Gravel Till Sand/Gravel Till Limestone Dolostone Evaporates/Shale/Limestone/ Dolostone Dolostone Dolostone Karstified dolostone Sandstone/Shale Shale

2.5  108–5.6  103 1  108–2.8  103 5  109–2.5  103 1  108–3.5  103 1.1  109–3.9  103 106–103 5  107–1  103 5  107–1  103 5  105 5  108/5  106 5  105 107 108

Please cite this article in press as: Hwang, H.-T., et al. A simple iterative method for estimating evapotranspiration with integrated surface/subsurface flow models. J. Hydrol. (2015), http://dx.doi.org/10.1016/j.jhydrol.2015.10.003

6

H.-T. Hwang et al. / Journal of Hydrology xxx (2015) xxx–xxx

Table 4 Manning’s roughness coefficients for the GRW land cover types (Engman, 1986; McCuen, 1989). Class name

Roughness coefficient (m1/3s)

Builtup (residential) Builtup (commercial/industrial) Row crops Small grains Forage Pasture/sparse forest Dense forest (deciduous) Dense forest (conifer) Dense forest (mixed) Plantation (mature) Open water Wetlands Extraction/bedrock/roads/beaches Golf courses Bare agricultural fields Other

0.015 0.012 0.1 0.2 0.13 0.07 0.6 0.6 0.6 0.1 0.03 0.05 0.011 0.45 0.3 0.15

Conservation Authority, 1999) and it contains 16 land cover classes. Among the 16 land cover classes, the areas for agriculture and forage cover 57% of total area and the urbanized areas occupy approximately 5% of the total area. Surface and subsurface flow parameters used in the HGS simulations are summarized in Tables 3 and 4. The 2D surface domain was discretized using various sizes of triangular elements, consisting of about 220,000 nodes and 437,000 elements. The subsurface domain consists of 14 element layers (13 subsurface material types) and was discretized using 3D triangular prism elements. The 3D finite element mesh consists of about 3.5  106 nodes and 6.6  106 elements. The average horizontal resolution of the 2D triangular mesh is about 200 m and it was refined along the surface drainage features at a resolution of 50 m or less (Fig. 5). Using the Climate Research Unit (CRU) climate datasets (Harris et al., 2014), an annual averaged total precipitation from 1979 to

1988 was applied as the rainfall boundary conditions. The annual average total precipitation ranges between 800 and 1100 mm/year in the watershed, while the annual average potential ET is estimated to range between 660 and 780 mm/year. A critical depth boundary condition is used around the entire surface domain. To perform the iterative method, it is important to classify the upstream and downstream relationships between gauging stations. The surface drainage network and the locations of the 16 stream flow gauging stations in the GRW are shown in Fig. 6. Note that gauging station 16 located at the outlet of the GRW is not an actual stream flow measurement site. The stream flow at station 16 used in this study is based on the ET simulation (i.e. reference simulation) results. Drainage areas associated with the gauging stations (Fig. 6a) were delineated based on the DEM data. As water balance in the upstream sub-watersheds can affect the downstream areas, it is important to define the relationships among the sub-watersheds. As shown in Fig. 6b, the stations 2, 3, 4, 11 and 15 are located along the Grand River, with stations 11 and 15 lying downstream of the Conestogo, Speed, and Nith River tributaries. Specifically, stream discharge at station 11 is directly affected by the discharge at stations 4, 7 and 10, which lie on the Grand, Conestogo, and Speed Rivers, respectively. Similarly, station 15 is located downstream of the Nith River sub-watershed (i.e. station 14) as well as the other Grand River tributaries. Surface water and variably-saturated subsurface flow was initially simulated in conjunction with ET processes to establish the benchmark case, which is compared to the actual ET estimated using the iterative method. To simulate the ET processes, the land cover classes were re-grouped as upland (dense forest and plantation areas), lowland (agricultural area such as row crops, small grains, forage, pasture, golf course, and bare agricultural fields), wetland (open water and wetland areas) and urban (built-up areas, roads, bedrock, and beaches). Table 5 lists the parameters used for the ET simulation. Using the given flow and ET parameters listed in Tables 3–5, the simulated steady-state stream flow rates and ET

Fig. 5. The GRW simulation domain and its finite element discretization.

Please cite this article in press as: Hwang, H.-T., et al. A simple iterative method for estimating evapotranspiration with integrated surface/subsurface flow models. J. Hydrol. (2015), http://dx.doi.org/10.1016/j.jhydrol.2015.10.003

H.-T. Hwang et al. / Journal of Hydrology xxx (2015) xxx–xxx

7

Fig. 6. Locations of stream flow gauge stations in the GRW (a) surface drainage networks; (b) subwatersheds. Note that station 16 is from the ET simulation (i.e. reference simulation) results.

Table 5 ET parameter values for the GRW simulation. Parameter

Upland

Lowland

Wetland

Urban

Evaporation depth (m) Root depth (m) LAI table (–)

1.0 2.9 5.12

1.0 1.5 4.5

1.0 0.5 0.0

1.0 0.1 0.4

Transpiration fitting parameters (–) C1 0.15 C2 0.1 C3 1.0

0.15 0.1 1.0

0.3 0.2 1.0

0.15 0.1 1.0

Transpiration limiting saturation (–) Wilting point 0.18 Field capacity 0.19 Oxic limit 1.0 Anoxic limit 1.0

0.18 0.19 1.0 1.0

0.18 0.19 1.0 1.0

0.18 0.19 1.0 1.0

Evaporation limiting saturation (–) Min 0.18 Max 0.19

0.18 0.19

0.18 0.19

0.18 0.19

Canopy storage (–) Initial interception (–)

0.0 0.0

0.0 0.0

0.0 0.0

0.0 0.0

values are obtained and used as reference for the model validations. It is to be noted that the performance of the reference model is not the main focus of this study as it primarily serves as a benchmark for which to compare the results from the simulations where ET was calculated using the newly proposed iterative method. The simulated actual ET distribution over the domain shows that the annual average ET rate can range from below 100 to up to 780 mm/year (maximum potential ET rate) and it is generally higher near streams and in the downstream areas where the water table intersects ground surface or is relatively shallow (Fig. 7). 4.2. Iterative sub-watershed evapotranspiration estimation To validate the accuracy and applicability of the iterative method presented in this study, the iterative method was applied

Fig. 7. Distribution of simulated average annual ET in the GRW. Note that station 16 is from the ET simulation (i.e. reference simulation) results.

to estimate the average ET rates in the 16 sub-watersheds within the GRW. The reference and iterative models are identical to each other, except that ET processes are only applied to the reference simulation. Because net precipitation cannot exceed effective precipitation, the initial net precipitation ratios (a) were assigned as 1.0 for all the sub-watersheds, which means that the net precipitation was assumed to be the same as the total precipitation. For convergence criteria, a 1.0% relative difference between simulated

Please cite this article in press as: Hwang, H.-T., et al. A simple iterative method for estimating evapotranspiration with integrated surface/subsurface flow models. J. Hydrol. (2015), http://dx.doi.org/10.1016/j.jhydrol.2015.10.003

8

H.-T. Hwang et al. / Journal of Hydrology xxx (2015) xxx–xxx Table 6 Spatially averaged net precipitation and actual evapotranspiration estimated from the iterative approach.

a

Sub-watershed

a

Estimated net precipitation (mm/year)

Estimated actual evapotranspiration (mm/year)

Dundalk Marsville Shand Dam Montrose Drayton Glen Allan St. Jacobs Amstrong Mills Guelph Cambridge Galt Nithburg Hamburg Canning Brantford GRW end

0.39 0.40 0.50 0.47 0.49 0.43 0.41 0.49 0.47 0.48 0.44 0.38 0.44 0.44 0.42 0.29

387.1 387.0 466.3 438.9 469.0 434.8 413.5 453.5 427.0 457.6 416.1 383.4 456.3 432.7 408.4 281.5

597.4 570.6 467.6 495.0 480.2 574.5 598.0 480.4 483.9 495.3 538.6 633.9 570.7 549.4 574.7 673.5

GRWa

0.41

390.6

573.8

Overall value for the GRW, which is taken into account the areal ratios.

Fig. 8. (a) Differences between synthetic and estimated surface flow rates at each iteration and (b) comparison between simulated and estimated flow rates after 3rd iteration.

and observed stream flow rates was assigned at each gauging station location; and during the simulations, convergence was achieved after the third iteration. The effectiveness of the iterative method is highlighted in Fig. 8, with Fig. 8a showing the flow rate difference in relation to the benchmark simulation at each iteration level for each of the 16 gauging station locations, and Fig. 8b showing a comparison of stream flow rates between the benchmark simulation and the iterative method, with the maximum difference being about 0.07 m3/s. Table 6 lists the spatially-averaged net precipitation and ET estimates generated with the iterative method, which are based on the a values obtained from the results of the 3rd iteration. The estimated ET for the sub-watersheds ranges from 467.6 to 673.5 mm/year (50.1% to 70.5% of the effective precipitation). Overall, the estimated annual net precipitation and ET rates are 390.6 and 573.8 mm/year, respectively. The spatially-averaged a over the domain is approximately 0.41, indicating that the net precipitation over the area is about 41% of the effective precipitation and the remainder (59%) of the precipitation considered lost to ET. Fig. 9 shows a spatially-averaged actual ET comparison between the benchmark simulation and

Fig. 9. Result comparison at each subwatershed between synthetic and estimated AET.

the iterative approach for the 16 sub-watersheds. Although the actual ET for the areas upstream of the Shand Dam and Brantford stations was slightly underestimated, there is a high level of agreement between the results from the benchmark simulations and the iterative approach simulations, with an overall R2 of 0.99. Although the iterative approach, combined with the integrated hydrologic model, could reproduce the observed stream flow at the gauging stations, a spatially uniform ET to precipitation ratio was assumed for each sub-watershed (see Fig. 10). When the results from the benchmark simulations (where the modelled ET is spatially variable) are compared to the results from the iterative approach simulations, it is found that the overall differences in both soil saturation and surface water depth are reasonably small (Fig. 11). Therefore, it can be surmised that at the sub-watershed scale within the GRW, the newly proposed iterative approach can effectively replicate the key elements of the simulations where

Please cite this article in press as: Hwang, H.-T., et al. A simple iterative method for estimating evapotranspiration with integrated surface/subsurface flow models. J. Hydrol. (2015), http://dx.doi.org/10.1016/j.jhydrol.2015.10.003

H.-T. Hwang et al. / Journal of Hydrology xxx (2015) xxx–xxx

9

Fig. 10. Distribution of estimated ET in the GRW. Note that station 16 is from the ET simulation (i.e. reference simulation) results.

ET was fully modelled by the integrated surface/subsurface flow model. It should be noted that this study suggests a simple and robust iterative approach to estimate the average rate of subwatershed scale ET from measured stream flow and is not intended to justify the use of average ET over large areas for all hydrologic modelling applications. However, in cases where the observational data required to calculate spatially variable ET are scarce, the newly proposed ET calculation method could prove valuable. 5. Summary and conclusions Evapotranspiration (ET) is one of the most important hydrologic processes in the water cycle for determining the balance of water between eco-hydrosphere and atmosphere. The rate of ET is difficult to measure even at a single location, not to mention its distribution over a watershed or sub-watershed scale. Alternatively, the average rate of ET over an area of interest can be estimated by analyzing the water budget when it is a single unknown and all the other components in the balance equation can be reasonably estimated, or measured. For example, when a watershed system can be approximated as a closed system with groundwater flow at steady-state, ET can be estimated relatively accurately by measuring precipitation and stream flow at all inlets and outlets. However, in a sub-drainage basin scale groundwater can flow in and out of the system, and thus quantifying all components of the water balance equation becomes more difficult. In these cases, the estimation of ET becomes more complex, for which tools such as numerical models can serve a valuable purpose. This study presents methodology for a combined use of the water balance equation and an integrated hydrologic model to estimate the average rates of ET for sub-watersheds associated with measured stream flows. A simple Newton iterative approach was developed to estimate the ratios of net precipitation to precipitation for sub-watersheds in order to minimize differences between measured and simulated stream flow rates. The method was first

Fig. 11. Distribution of absolute difference between synthetic and estimated (a) saturation and (b) surface water depth.

applied to a synthetic watershed for illustration purposes and for proof-of-concept demonstration. This proof-of-concept example showed that the convergence for the iterative method could be achieved with a reasonable number of iterations (3–5 depending on the initial guess). It was also shown that the simplifying assumptions were not limiting factors for the simple virtual watershed used for the proof of concept.

Please cite this article in press as: Hwang, H.-T., et al. A simple iterative method for estimating evapotranspiration with integrated surface/subsurface flow models. J. Hydrol. (2015), http://dx.doi.org/10.1016/j.jhydrol.2015.10.003

10

H.-T. Hwang et al. / Journal of Hydrology xxx (2015) xxx–xxx

The applicability of the method was then tested for a highlycharacterized real-world watershed located in Southern Ontario, Canada (the 7000 km2 Grand River Watershed). The surface drainage network for the watershed consists of one main stream and its multiple tributaries, and the average stream flow rates measured at 16 gauging stations were utilized to estimate the ratios of ET to precipitation in 16 sub-watersheds. As the gauge stations are located along the major river and also at tributaries, the relation among the sub-watersheds associated with the stream gauging stations is a multi-path chain. Only three iterations were required for convergence that resulted in less than 1% relative error in stream flow calibration at all stations. Comparison between the calibrated surface/subsurface flow model that explicitly simulated ET and the calibrated model that utilized the newly proposed iterative method shows that the iterative approach (without explicitly simulating ET) could effectively approximate ET for the system using 16 spatially distributed stream flow measurements. Although this approach is applicable to any well-characterized steady-state hydrologic system where an integrated surface/ subsurface flow model is available, it should be noted that this approach can only provide spatially-averaged estimates of static ET when the model is calibrated to measured stream flows, in order to satisfy the requirements of a water budget based analysis. The validity of the calibrated model will thus depend on the density and quality of the data used and also the validity of the integrated hydrologic model. Because temporal changes in ET can significantly impact short and long term water budgets in the hydrosphere, and moisture distribution in atmosphere, extending the iterative methodology presented here to provide estimations of temporally-variable ET under transient conditions will be the next step in this work. Acknowledgements The authors gratefully thank Grand River Conservation Authority for providing the datasets for the Grand River Watershed. The simulations discussed in this paper were performed on the General Purpose Cluster at the SciNet High Performance Computing facility of the University of Toronto. Funding was partially provided by grants from an IBM Southern Ontario Smart Computing Innovation Platform (SOSCIP) Post-doctoral Fellowship; and from the Agricultural Research Institute of Ontario and the Ontario Ministry of Agriculture, Food and Rural Affairs through the New Directions Research Program. References Ala-aho, P., Rossi, P.M., Isokangas, E., Kløve, B., 2015. Fully integrated surface– subsurface flow modelling of groundwater–lake interaction in an esker aquifer: model verification with stable isotopes and airborne thermal imaging. J. Hydrol. 522, 391–406. Aquanty Inc., 2015. HydroGeoSphere. A Three-dimensional Numerical Model Describing Fully-integrated Subsurface and Surface Flow and Solute Transport. Waterloo, ON. AquaResource Inc., 2009. Integrated Water Budget Report, Grand River Watershed. Submitted to the Grand River Conservation Authority, Cambridge, ON, Canada. Bearup, L.A., Maxwell, R.M., Clow, D.W., McCray, J.E., 2014. Hydrological effects of forest transpiration loss in bark beetle-impacted watersheds. Nat. Clim. Change 4 (6), 481–486. Bolger, B.L., Park, Y.-J., Unger, A.J.A., Sudicky, E.A., 2011. Simulating the predevelopment hydrologic conditions in the San Joaquin Valley, California. J. Hydrol. 411 (3–4), 322–330. Buerge, I.J., Buser, H.-R., Kahle, M., Müller, M.D., Poiger, T., 2009. Ubiquitous occurrence of the artificial sweetener acesulfame in the aquatic environment: an ideal chemical marker of domestic wastewater in groundwater. Environ. Sci. Technol. 43 (12), 4381–4385. Condon, L.E., Maxwell, R.M., 2014. Feedbacks between managed irrigation and water availability: diagnosing temporal and spatial patterns using an integrated hydrologic model. Water Resour. Res. 50 (3), 2600–2616. Dastorani, M., Poormohammadi, S., 2012. Evaluation of water balance in a mountainous upland catchment using SEBAL approach. Water Resour. Manage. 26 (7), 2069–2080.

Davis, S.L., Dukes, M.D., 2010. Irrigation scheduling performance by evapotranspiration-based controllers. Agric. Water Manage. 98 (1), 19–28. De Schepper, G., Therrien, R., Refsgaard, J.C., Hansen, A.L., 2015. Simulating coupled surface and subsurface water flow in a tile-drained agricultural catchment. J. Hydrol. 521, 374–388. Delfs, J.-O., Wang, W., Kalbacher, T., Singh, A., Kolditz, O., 2013. A coupled surface/subsurface flow model accounting for air entrapment and air pressure counterflow. Environ. Earth Sci. 69 (2), 395–414. Ebel, B.A., Loague, K., 2006. Physics-based hydrologic-response simulation: seeing through the fog of equifinality. Hydrol. Process. 20 (13), 2887–2900. Engman, E., 1986. Roughness coefficients for routing surface runoff. J. Irrig. Drain. Eng. 112 (1), 39–53. Frey, S.K., Topp, E., Ball, B.R., Edwards, M., Gottschall, N., Sunohara, M., Zoski, E., Lapen, D. R., 2013. Tile drainage management influences on surface-water and groundwater quality following liquid manure application. J. Environ. Qual. 42 (3), 881–892. Goderniaux, P., Brouyère, S., Fowler, H.J., Blenkinsop, S., Therrien, R., Orban, P., Dassargues, A., 2009. Large scale surface–subsurface hydrological model to assess climate change impacts on groundwater reserves. J. Hydrol. 373 (1–2), 122–138. Gowda, P., Chavez, J., Colaizzi, P., Evett, S., Howell, T., Tolk, J., 2008. ET mapping for agricultural water management: present status and challenges. Irrig. Sci. 26 (3), 223–237. Grand River Conservation Authority, 1999. Land Cover Data. Grand River Conservation Authority, Waterloo, ON, . Hanson, R.L., 2009. Evapotranspiration and droughts. In: Paulson, R.W., Chase, E.B., Roberts, R.S., Moody, D.W. (Eds.), Compilers, National Water Summary 1988– 89-hydrologic Events and Floods and Droughts: U.S. Geological Survey WaterSupply Paper 2375, pp 99–104. Harris, I., Jones, P.D., Osborn, T.J., Lister, D.H., 2014. Updated high-resolution grids of monthly climatic observations – the CRU TS3.10 Dataset. Int. J. Climatol. 34 (3), 623–642. Hassan, S.M.T., Lubczynski, M.W., Niswonger, R.G., Su, Z., 2014. Surface– groundwater interactions in hard rocks in Sardon Catchment of western Spain: an integrated modeling approach. J. Hydrol. 517, 390–410. Hornberger, G.M, 1998. Elements of Physical Hydrology. JHU Press. Hwang, H.T., Park, Y.J., Sudicky, E.A., Forsyth, P.A., 2014. A parallel computational framework to solve flow and transport in integrated surface–subsurface hydrologic systems. Environ. Model. Softw. 61, 39–58. Kollet, S.J., Maxwell, R.M., 2006. Integrated surface–groundwater flow modeling: a free-surface overland flow boundary condition in a parallel groundwater flow model. Adv. Water Resour. 29 (7), 945–958. Kristensen, K.J., Jensen, S.E., 1975. A model for estimating actual evapotranspiration from potential evapotranspiration. Nord. Hydrol. 6 (3), 170–188. Kundzewicz, Z.W., Mata, L.J., Arnell, N.W., Döll, P., Kabat, P., Jiménez, B., Miller, K.A., Oki, T., Sen, Z., Shiklomanov, I.A., 2007. Freshwater resources and their management. Climate change 2007: impacts, adaptation and vulnerability. In: Parry, M.L., Canziani, O.F., Palutikof, J.P., van der Linden, P.J., Hanson, C.E. (Eds.), Contribution of Working Group II to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge University Press, Cambridge, UK, pp. 173–210. Lake Erie Source Protection Region, 2008. Grand River Watershed Characterization Report. Grand River Conservation Authority, Waterloo, ON. Li, L., Lambert, M.F., Maier, H.R., Partington, D., Simmons, C.T., 2015. Assessment of the internal dynamics of the Australian Water Balance Model under different calibration regimes. Environ. Model. Softw. 66, 57–68. Li, Q., Unger, A.J.A., Sudicky, E.A., Kassenaar, D., Wexler, E.J., Shikaze, S., 2008. Simulating the multi-seasonal response of a large-scale watershed with a 3D physically-based hydrologic model. J. Hydrol. 357 (3–4), 317–336. Maneta, M.P., Schnabel, S., Wallender, W.W., Panday, S., Jetten, V., 2008. Calibration of an evapotranspiration model to simulate soil water dynamics in a semiarid rangeland. Hydrol. Process. 22 (24), 4655–4669. Maxwell, R.M., Chow, F.K., Kollet, S.J., 2007. The groundwater–land-surface– atmosphere connection: soil moisture effects on the atmospheric boundary layer in fully-coupled simulations. Adv. Water Resour. 30 (12), 2447–2466. McCuen, R.H., 1989. Hydrologic Analysis and Design. Prentice-Hall Englewood Cliffs, NJ. Niu, G.-Y., Paniconi, C., Troch, P.A., Scott, R.L., Durcik, M., Zeng, X., Huxman, T., Goodrich, D.C., 2014. An integrated modelling framework of catchment-scale ecohydrological processes: 1. Model description and tests over an energylimited watershed. Ecohydrology 7 (2), 427–439. Panday, S., Huyakorn, P.S., 2004. A fully coupled physically-based spatiallydistributed model for evaluating surface/subsurface flow. Adv. Water Resour. 27 (4), 361–382. Partington, D., Brunner, P., Simmons, C.T., Therrien, R., Werner, A.D., Dandy, G.C., Maier, H.R., 2011. A hydraulic mixing-cell method to quantify the groundwater component of streamflow within spatially distributed fully integrated surface water–groundwater flow models. Environ. Model. Softw. 26 (7), 886–898. Pérez, A.J., Abrahão, R., Causapé, J., Cirpka, O.A., Bürger, C.M., 2011. Simulating the transition of a semi-arid rainfed catchment towards irrigation agriculture. J. Hydrol. 409 (3–4), 663–681. Quinn, P.F., Beven, K.J., 1993. Spatial and temporal predictions of soil moisture dynamics, runoff, variable source areas and evapotranspiration for plynlimon, mid-wales. Hydrol. Process. 7 (4), 425–448. Rassam, D.W., Peeters, L., Pickett, T., Jolly, I., Holz, L., 2013. Accounting for surface– groundwater interactions and their uncertainty in river and groundwater models: a case study in the Namoi River, Australia. Environ. Model. Softw. 50, 108–119.

Please cite this article in press as: Hwang, H.-T., et al. A simple iterative method for estimating evapotranspiration with integrated surface/subsurface flow models. J. Hydrol. (2015), http://dx.doi.org/10.1016/j.jhydrol.2015.10.003

H.-T. Hwang et al. / Journal of Hydrology xxx (2015) xxx–xxx Robertson, W.D., Van Stempvoort, D.R., Solomon, D.K., Homewood, J., Brown, S.J., Spoelstra, J., Schiff, S.L., 2013. Persistence of artificial sweeteners in a 15-yearold septic system plume. J. Hydrol. 477, 43–54. Rodell, M., Famiglietti, J.S., 2002. The potential for satellite-based monitoring of groundwater storage changes using GRACE: the High Plains aquifer, Central US. J. Hydrol. 263 (1–4), 245–256. Rodell, M., Famiglietti, J.S., Chen, J., Seneviratne, S.I., Viterbo, P., Holl, S., Wilson, C.R., 2004. Basin scale estimates of evapotranspiration using GRACE and other observations. Geophys. Res. Lett. 31 (20), L20504. Sellers, P.J., Randall, D.A., Collatz, G.J., Berry, J.A., Field, C.B., Dazlich, D.A., Zhang, C., Collelo, G.D., Bounoua, L., 1996. A revised land surface parameterization (SiB2) for atmospheric GCMS. Part I: model formulation. J. Clim. 9 (4), 676–705.

11

Van Stempvoort, D.R., Robertson, W.D., Brown, S.J., 2011. Artificial sweeteners in a large septic plume. Ground Water Monit. Rem. 31 (4), 95–102. VanderKwaak, J.E., Loague, K., 2001. Hydrologic-response simulations for the R-5 catchment with a comprehensive physics-based model. Water Resour. Res. 37 (4), 999–1013. Xue, B.-L., Wang, L., Li, X., Yang, K., Chen, D., Sun, L., 2013. Evaluation of evapotranspiration estimates for two river basins on the Tibetan Plateau by a water balance method. J. Hydrol. 492, 290–297. Yi, Y., Gibson, J.J., Hélie, J.-F., Dick, T.A., 2010. Synoptic and time-series stable isotope surveys of the Mackenzie River from Great Slave Lake to the Arctic Ocean, 2003 to 2006. J. Hydrol. 383 (3–4), 223–232.

Please cite this article in press as: Hwang, H.-T., et al. A simple iterative method for estimating evapotranspiration with integrated surface/subsurface flow models. J. Hydrol. (2015), http://dx.doi.org/10.1016/j.jhydrol.2015.10.003