e c o l o g i c a l m o d e l l i n g 2 0 4 ( 2 0 0 7 ) 503–515
available at www.sciencedirect.com
journal homepage: www.elsevier.com/locate/ecolmodel
A water balance model for rainfed lowland rice fields emphasising lateral water movement within a toposequence M. Tsubo a,∗ , S. Fukai a , T.P. Tuong b , M. Ouk c a b c
School of Land, Crop and Food Sciences, The University of Queensland, Brisbane, Qld 4072, Australia Crop and Environmental Sciences Division, International Rice Research Institute, DA PO Box 7777, Metro Manila, Philippines Cambodian Agricultural Research and Development Institute, Phnom Penh, Cambodia
a r t i c l e
i n f o
a b s t r a c t
Article history:
In the Mekong region, most rice paddies lie in gently sloping lowlands. An accurate esti-
Received 8 March 2006
mate of paddy water availability is crucial for modelling rice (Oryza sativa L.) productivity
Received in revised form
in such environments. However, modelling of water balance in the sloping lowlands can
22 January 2007
be difficult due to problems in estimating lateral water movement from high to low posi-
Accepted 1 February 2007
tions in the toposequence. This paper describes a semi-empirical model for estimating net
Published on line 28 February 2007
lateral water flow along a toposequence of rice fields. Net lateral water flow is separated into three sub-components: (i) lateral seepage through the bunds from field to field; (ii) sur-
Keywords:
face runoff over the bunds from field to field; (iii) water run-on from the catchment above
Modelling
the toposequence. The lateral seepage is estimated using the Dupuit equation for steady
Oryza sativa
unconfined flow. Surface runoff over the bund is calculated as excess water depth above
Percolation
the bund height, while run-on from the catchment is calculated using a rainfall–runoff rela-
Runoff
tionship. The other water balance components used in the model are rainfall (measured),
Seepage
paddy evapotranspiration (ET), and downward water flow (within the field and under the bund). Paddy ET is estimated using the FAO crop ET model. The downward water flow is empirically determined for soil texture groups. Simulation results showed that field water levels computed by the model varied between the bund height (maximum) and the hardpan level (minimum). The daily change in the water level was satisfactorily validated using data collected from a field experiment in a toposequence in Cambodia. The component governing the water balance in the higher toposequence positions was the downward water flow (45% of seasonal rainfall at the top of the toposequence). However, the water balance in the lower positions was affected largely by the net lateral water flow (21% of seasonal rainfall at the bottom of the toposequence). Seasonal ET was 67–69% of seasonal rainfall. © 2007 Elsevier B.V. All rights reserved.
1.
Introduction
Drought is a major problem in the rainfed lowland rice (Oryza sativa L.) growing areas of Southeast Asia, particularly in
Cambodia, Laos and Thailand (the Mekong region). Over 70% of the rice area in the Mekong region is classified as belonging to the rainfed lowland rice ecosystem; of this, more than 70% belongs to the drought-prone and drought-and-
∗ Corresponding author. Present address: 1390 Hamasaka, Arid Land Research Center, Tottori University, Tottori 680 0001, Japan. Tel.: +81 857 21 7037; fax: +81 857 21 7037. E-mail address:
[email protected] (M. Tsubo). 0304-3800/$ – see front matter © 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.ecolmodel.2007.02.001
504
e c o l o g i c a l m o d e l l i n g 2 0 4 ( 2 0 0 7 ) 503–515
Nomenclature catchment area (m2 ); ‘above’ denotes AC above toposequence paddy field area (m2 ) AP C topsoil clay content (%) D hardpan depth (m) Dtotal total downward water flow (m day−1 ) downward water flow within field (m day−1 ); Dunder i.e., percolation through hardpan Dwithin downward water flow under bund (m day−1 ); i.e., under-bund percolation loss ETpaddy paddy evapotranspiration (m day−1 ); ‘stress’ denotes ETpaddy under water deficit reference evapotranspiration (m day−1 ) ETo fx fractional distance of toposequence; ‘centre’ denotes fx at the centre of field fractional width of toposequence; ‘centre’ fy denotes fy at the centre of field fz fractional elevation of toposequence; ‘centre’ denotes fy at the centre of field H bund height (m) water level of earth dam reservoir; ‘high’ and Hdam ‘low’ denote high and low H, respectively HPW water head from hardpan (m) K hydraulic conductivity of bund (m day−1 ); ‘dam’ denotes k for earth dam crop coefficient Kc Ks plant water availability coefficient; ‘d’ denotes Ks when free water surface is at hardpan Lin lateral water inflow (m day−1 ); ‘through’ and ‘over’ denote Lin through and over bund, respectively Lnet net lateral water flow (m day−1 ) lateral water outflow (m day−1 ); ‘through’ and Lout ‘over’ denote Lout through and over bund, respectively Lrunon run-on of water from catchment (m day−1 ) N number of paddy fields along toposequence water flow rate through earth dam per unit qdam bund length (m3 day−1 ) qthrough lateral seepage flow rate per unit bund length (m3 day−1 ) under-bund percolation rate per unit bund qunder length (m3 day−1 ) R rainfall (m day−1 ) minimum rainfall required for runoff to occur Rmin (m day−1 ) S slope of the saturation vapour pressure curve S change in water depth (m day−1 ) w bund thickness (m); ‘dam’ denotes w for earth dam Winput input water depth to field (m) WG groundwater table depth (m) Woutput output water depth to field (m) field water level above hardpan (m) WP xcatchment length of catchment area above toposequence (m) AC
X ytop Y Z
total distance of toposequence (m) width of toposequence at the top (m) width of toposequence at the bottom (m) elevation of toposequence at the top relative to the bottom (m)
Greek letters ˛ runoff coefficient for soil infiltration ˇ constant for determining minimum rainfall required for runoff to occur (mm) ı coefficient for the Priestley and Taylor (1972) evapotranspiration model ε coefficient for reducing downward water flow porosity of soil (m3 m−3 ) coefficient for the slope of toposequence
submergence-prone environments (Mackill et al., 1996; IRRI, 2005). Therefore, approximately half of the rice growing area of the Mekong region potentially encounters drought. Even though total rainfall might be adequate for rice cultivation, poor distribution during the growing season can sometimes lead to crop failure. Simulation models have been employed to analyse crop responses to water availability in rice cropping systems, and the model outputs, including growth and yield, are reasonably acceptable. The water environments characterized at a single field level, provide values for the water balance components such as deep percolation and lateral water movement (Fukai et al., 1995, 2000). However, none of the simulation models available adequately demonstrate spatial variability of the crop responses in water-limiting environments within a watershed (a few kilometres in distance), though there are some models dealing with large catchments (Becu et al., 2003). “Toposequence” is defined as a sequence of paddies lying on sloping land, normally less than 1 km long. The upper and lower paddies in a toposequence tend to be classified as providing drought- and submergence-prone environments, respectively (Wade et al., 1999), reflecting an association between toposequence position and the nature of the rice environment. Several water balance models for toposequences have been developed. For example, PADIWATER (Bolton and Zandstra, 1981) is a simple water balance model dealing with surface runoff along a toposequence, in which surface runoff from one paddy field is treated as water input to the next paddy field below. More complex models for water movement in paddy fields have been developed on the basis of Darcy’s law (McMennamy and O’Toole, 1983; Odhiambo and Murty, 1996; Bouman et al., 2001; Chen et al., 2002). The one-dimensional models, such as ORYZA2000 (Bouman et al., 2001), CERES-Rice (Singh et al., 1993; Mohmood, 1998) estimates only deep percolation, whereas for the sloping landscapes, two- and three-dimensional models such as FEMWATER are employed to estimate both the vertical and horizontal water flows (Chen et al., 2002). These models require detail information on the hydraulic conductivity of paddy soils; such information is usually not readily obtainable. Therefore, semiempirical models may be required for estimates of water
505
e c o l o g i c a l m o d e l l i n g 2 0 4 ( 2 0 0 7 ) 503–515
balance in a toposequence. In model development, water balance components that represent lateral water movement in the sloped lowland, particularly lateral seepage from field to field, should be considered, and the Dupuit equation for steady unconfined flow (an analytical solution for a two-dimensional flow domain) can be applied for the lateral water flow problem. The timing of the disappearance of standing water in paddies differs with toposequence positions. This toposequential effect might explain why rice productivity varies across the landscapes (Tsubo et al., 2006). In current drought-tolerant rice breeding programmes in the Mekong region, the importance of quantifying water availability along a toposequence, is often not realized. Part of the reason for this is on account of the complexity of lateral water movement in the sloping lowlands, as it is not easy to measure water loss by lateral movement. An earlier study (Tsubo et al., 2006) has shown that lateral water flows in and out of paddies throughout the growing season, and this lateral water inflow–outflow processes can be quantified in relation to rainfall. However, in addition to direct rainwater catchment, water runoff from adjoining nonpaddy areas (or catchment) can be an additional water input to a paddy field. Therefore, quantifying lateral water movement associated with rainfall is not straightforward. In modelling water balance in sloping lowlands, the impact of water runon from adjoining catchment areas on water balance needs to be taken into account. Suzuki et al. (2003) recently formulated a two-layered “tank” model for perched water and groundwater balances along a toposequence in Northeast Thailand, in which each component of the water balance was parameterized. This model deals with run-on of water from the catchment above the uppermost field in the paddy toposequence; however, it does not take account of possible inflow from catchments to the side of the toposequence. Furthermore, when perched water is present in a paddy field, under-bund percolation should be taken into consideration in the water balance, as it can account for a significant amount of the water lost from the field (Tuong et al., 1994). The objective of this study was, therefore, to develop a model for water balance along a paddy toposequence, with particular reference to lateral water movement from the catchment to higher paddy fields and then to lower paddy fields, applying a geometrical approach to take into consideration the differing dimensions of the individual components of the toposequence. The model developed was based on data collected from a field experiment in Kampong Cham province in Cambodia (12◦ 00 N, 105◦ 27 E), and was then subjected to a sensitivity analysis.
toposequence was about 0.4 km; the slope of the toposequence ranged from 1% at lower positions, to 3% in the highest positions. Two different types of free water level measurements were taken: a short PVC tube was used to measure the water level above the hardpan (the perched water table), and a long PVC tube for groundwater level. The water levels were measured on a daily basis at five toposequence positions (top, middle-1, middle-2, middle-3 and bottom). The apparatus used to measure water levels is described in detail elsewhere (Tsubo et al., 2005). Daily rainfall and pan evaporation were also recorded at the site.
2.2.
Model description
The watershed on which the model is based included the catchment generating the water runoff to the toposequence. The specifications of field and catchment areas were determined for each toposequence position, to quantify lateral water movement from field-to-field, and from catchment-tofield. For each toposequence position, the soil surface was defined as reference elevation (i.e., zero) for the purpose of the computations. The bund height, h (a positive value), and the hardpan depth, d (a negative value), were assumed to be constant relative to the soil surface, for each toposequence position. Other important assumptions are listed in Table 1. A water balance equation is used to estimate free water level above the hardpan (WP ) along the toposequence. When perched water is present, the groundwater level relative to the soil surface (WG ), which appears below the hardpan (WG < d), is differentiated from WP . When WP disappears from the topsoil (the soil layer above the hardpan), only WG is present in the subsoil (the soil below the hardpan). A description of the water balance equation employed is as follows. The change in the water depth at each toposequence position, S, is written as: S = R − ETpaddy − Dtotal − Lnet ,
mm of water
(1)
where R is rainfall, ETpaddy the evapotranspiration from the paddy field, Dtotal the total downward water flow, and Lnet is the net lateral water flow. Positive and negative S indicate water gain and loss, respectively. The spatial distribution of rainfall over the toposequence is assumed to be uniform, while the other components (i.e., ETpaddy , Dtotal and Lnet ) are calculated for each toposequence position.
2.2.1. Expression of paddy field positions in a toposequence
2.
Materials and methods
2.1.
Field experiment
The field experiment was carried out in a farmers’ toposequence of paddy fields, which were cultivated to several rice varieties in the 2004 wet-season (August–October). The soil in the locality was a loam with a 26% clay content. The paddy fields lay on gently sloping land, with hardpans at a depth of 20–30 cm and 10–40 cm high bunds around the fields. The horizontal distance from the top to the bottom of the
The size and shape of toposequences vary between locations. In order to standardise the toposequences for the purposes of the study, the distance, width and elevation of each toposequence (x, y and z, respectively) are expressed in fractional coordinates: i.e., fractional distance (fx = x/X where X is the projected distance of a toposequence on the xy-plane), fractional width (fy = y/Y where Y is the width of the toposequence at the bottom) and fractional elevation (fz = z/Z where Z is the elevation of the toposequece at the top, relative to the bottom). A typical toposequence can be depicted as has having an increasing slope and decreasing width, with increased
506
e c o l o g i c a l m o d e l l i n g 2 0 4 ( 2 0 0 7 ) 503–515
Table 1 – Assumptions made and inputs used for the model test for the study site (a) Assumptions Bund height and hardpan depth are constant relative to the soil surface along the toposequence The spatial distribution of rainfall over the toposequence is uniform An individual toposequence is made up of a number of level paddy fields on sloping land, and the paddy fields are parallel rectangular planes and horizontally separated by an equal distance The slope of catchment area above paddy fields dose not affect the coefficients ˛ (ratio of runoff to rainfall after the onset of run-off) and ˇ (a constant characterizing the minimum amount of rainfall required for runoff occurs) Hardpan is impermeable (to satisfy Dupuit-Forchheimer assumptions) When free water level is below the soil surface, the downward water flow declines linearly with the drop of the water level from the soil surface to the level of the hardpan The depth of paddy standing water does not affect downward water movement When free water level is below the soil surface, paddy evapotranspiration declines linearly with the drop of the water level from the soil surface to the level of the hardpan (b) Inputs Topsoil clay content, C Hardpan depth, D Bund height, H Hydraulic conductivity of bund, k Plant water availability coefficient for WP = d, Ks,d Number of paddy fields along toposequence, n Under-bund percolation rate per unit bund length, qunder Bund thickness, w Length of catchment area above toposequence, xcatchment Total distance of toposequence, X Width of toposequence at the top, ytop Width of toposequence at the bottom, Y Elevation of toposequence from the bottom, Z Runoff coefficient for soil infiltration, ˛ Constant for minimum rainfall required for runoff, ˇ Coefficient for the slope of toposequence, Porosity of paddy soil,
26% −0.25 m 0.25 m 0.8 m day−1 (for a loamy soil) 0.5 17 0.14 m3 day−1 (for a loamy soil) 0.5 m 50 m 372 m 45 m 124 m 4.6 m 0.3 3 mm 0.78 0.45 m3 m−3 (for a loamy soil)
elevation; the model described here has been developed for this type of toposequence (Fig. 1a). Given fz = 1 at fx = 0 for the top of the toposequence, and fz = 0 at fx = 1 for the bottom on the xz-plane, the following empirical function can be used to describe the change in slope with increased elevation: fz =
(1 − fx ) + fx
(2)
where is a coefficient for the slope of the toposequence ( > 0). An example of the use of Eq. (2) is illustrated in Fig. 2a ( = 0.78, r2 = 0.98 for one experimental site in Kampong Cham). If the slope is constant, fz decreases linearly with increases in fx (i.e., fz = 1 − fx ). By defining ytop as the width of the toposequence at the top (i.e., fy = ytop /Y at fx = 0 and fy = 1 at fx = 1 on the xy-plane), the linear relationship between fractional distance and width of the toposequence
Fig. 1 – Model toposequence in xyz co-ordinations (increased slope % and decreased width with increase in elevation): (a) continuous and (b) discrete.
507
e c o l o g i c a l m o d e l l i n g 2 0 4 ( 2 0 0 7 ) 503–515
Fig. 2 – Fractional distance, width and elevation of the toposequence in Kampong Cham, Cambodia ((a) xz-plane and (b) xy-plane), and positions of free water level (WP ) measurements along the toposequence ((c).
is given as: fy =
ytop ytop + 1− Y Y
fx
(3)
The ratio of ytop to Y (ytop /Y) is 0.36 for the Kampong Cham site (Fig. 2b). In addition, for a constant width, with increases in elevation, fy remains unchanged with increases in fx (i.e., fy = 1 for ytop = Y). An individual toposequence is assumed to be made up of a number of level paddy fields on sloping land. The soil surface of the paddy fields is considered to comprise parallel rectangular planes perpendicular to the x-axis and located at the xy-plane (see Fig. 1b for 10 paddy fields). Assuming that all the paddy fields are horizontally separated by an equal distance (17 toposequence positions for the Kampong Cham case; Fig. 2c) as the field size is similar between higher and lower positions in the toposequence concerned, the fractional distance at the centre of the ith position from the top of the toposequence, fx,centre(i) , is defined as: fx,centre(i) =
1 n
i−
1 2
(4)
where n is the total number of paddy fields along the toposequence (i = 1 for the top position, n for the bottom position). The fractional elevation and width at the centre of the ith position (fz,centre(i) and fy,centre(i) respectively) are therefore determined by substituting Eq. (4) into Eqs. (2) and (3), respectively. A geometrical relationship on the xy-plane is also used to determine field and catchment areas for a given rectangular area of the watershed concerned (Fig. 3). The field area at the ith position, AP(i) , and the catchment area for the ith paddy field, AC(i) , are given by: AP(i) =
X (Yfy,centre(i) ) n
(5)
and AC(i) = (xcatchment + Xfx,centre(i) ) ×
Y − ytop n
(6)
where xcatchment is the length of the catchment above the toposequence. Also, as shown in Fig. 3, there is an additional catchment area for the top paddy field in the toposequence (i.e., AC,above = xcatchment × ytop for n = 1).
508
e c o l o g i c a l m o d e l l i n g 2 0 4 ( 2 0 0 7 ) 503–515
Fig. 4 – Graphical expression of the relationship between runoff and rainfall in catchment areas (˛ and ˇ: runoff coefficients; ˇ/˛: the minimum rainfall amount required for the run-on to occur).
ratio of catchment to paddy field area as follows: Lrunon =
Fig. 3 – Graphical expression of the relationship between paddy and catchment areas (AP and AC , respectively) on the xy-plane.
2.2.2.
Net lateral water flow (Lnet )
The net lateral water flow (Lnet, positive indicates net outflow, negative net inflow) involves the inflow and outflow processes in the paddy field (Lin and Lout , respectively), and the run-on of water from the catchment to the paddy field (Lrunon ): Lnet = (Lout,over − Lin,over ) + (Lout,through − Lin,through ) − Lrunon (7) where ‘through’ and ‘over’ denote the lateral water movement through and over the bund, respectively. The lateral water outflow (Lout,through ) includes ponding water flow into the bund adjacent to the next field below, and soil water movement taking place in the cross-section between the hardpan level and the soil surface.
2.2.2.1. Lateral water flow from catchment to paddy field (Lrunon ). Run-on of water occurs only when a paddy field is directly adjacent to a catchment (Fig. 3). In general, runoff in a catchment occurs when rainfall exceeds infiltration into the soil. Run-on of water to a given paddy area from its catchment depends on the catchment size. Therefore, Lin,runon can be expressed as a proportion of rainfall by multiplying it by the
AC (˛R − ˇ) AP
(8)
where ˛ is the ratio of runoff to rainfall after the onset of run-off (˛ > 0), and ˇ (mm) is a constant for determining the minimum amount of rainfall required for runoff to occur (ˇ ≥ 0). The minimum rainfall (Rmin ) is determined from ˛ and ˇ (i.e., Rmin = ˇ/˛; Fig. 4). When a catchment receives less rainfall than Rmin (assumed to be 10 mm in the present study), there is no runoff from the catchment (i.e., Lin,runon = 0 for R ≤ ˇ/˛). The runoff–rainfall relationship is dependent on precipitation, soil and vegetation types. Suzuki et al. (2003) cites ˛ = 0.3 for a toposequence on gentle sloping land in Northeast Thailand; the lowland environment where Suzuki and colleagues did their water balance studies, had similar precipitation, soil and vegetation types to the site in the present study. The slope of the catchment area affects the coefficients ˛ and ˇ, but for the Kampong Cham case, Cambodia, this effect was probably minimal as the slope was very gentle (less than 3%).
2.2.2.2. Lateral water flow through bunds between paddy fields (Lout,through and Lin,through ). The Dupuit equation for steady unconfined flow can be applied for water flow problems in a two-dimensional domain. The assumptions (DupuitForchheimer assumptions) are: (a) there are no hydraulic gradients in the vertical dimension; (b) the hydraulic gradient in the horizontal dimension equals the slope of the water table (Reddi, 2003). These assumptions are valid for mild sloping land and shallow water depth. In the present study, a water balance model is being developed for such conditions (less than 5% slope, maximum field standing water level that is equal to the bund height). An application of the Dupuit equation is found in modelling water flow through earth dam (Fig. 5a). The Dupuit equation for estimating water flow rate through the dam (qdam ) is given by: qdam =
kdam 2 (H2 − Hdam,low ) 2wdam dam,high
(9a)
509
e c o l o g i c a l m o d e l l i n g 2 0 4 ( 2 0 0 7 ) 503–515
head above the hardpan (i.e., HPW = WP − d ≤ h − d). k values are very dependent on soil texture. For example, Tuong et al. (1994) reported 0.5 m day−1 for a clayey soil, and Walker and Rushton (1984) gave an estimate of 0.8 m day−1 for a loamy soil. For sandy soils, k values can be assumed to be greater than 1 m day−1 . When HPW(i+1) is positioned below the hardpan level of the ith position, there may be no effect of HPW(i+1) on qthrough(i) [i.e., HPW(i+1) − Z(fz(i) − fz(i+1) ) = 0 for HPW(i+1) ≤ Z(fz(i) − fz(i+1) ).]. For the bottom position in the toposequence (i = n), q(n) is an empirically given value as the under-bund percolation rate per unit length of the bund (qunder ), as explained below. In the case where the standing water disappears, Lout,through is very small, and is assumed to be nil (i.e., Lout,through = 0 for WP ≤ 0). Lout,through at the ith position is calculated as: Lout,through(i) =
Yfy centre(i) AP(i)
qthrough(i)
(10)
where Yfy,centre is the bund length parallel to the y-axis. Lout,through is regarded as input water to the paddy field immediately below, so Lin,through is estimated by weighting Lout,through by the area ratio of a paddy field to that of the next paddy field below (for i > 1):
Fig. 5 – Cross-section of field-bund-field showing the relationship of water head (HPW ) between paddies: (a) earth dam and (b) rice paddy bund. w: bund thickness; h: bund height; d: hardpan depth; Z(fz(i) − fz(i+1) ): elevation difference between paddies. Shaded areas indicate soil.
where kdam is the horizontal hydraulic conductivity of the dam, wdam the dam width, and Hdam is the water level of reservoirs (‘high’ and ‘low’ denote high and low water head, respectively). In the toposequence concerned, the Dupuit-Forchheimer assumptions are invalid unless the subsoil of the paddy field becomes fully saturated (i.e., WG ≥ d) (Walker and Rushton, 1984). In the Kampong Cham study, it was assumed that the hardpan was impermeable, as the downward water flow through the hardpan (less than 2 mm day−1 ) was very small relative to lateral water movement into the bund. When downward water flow under the bund (Dunder ), which is explained below (Section 2.2.3.2), occurs from the (i + 1)th paddy field to the bund adjacent to the ith paddy field, or when no perched water is developed in the (i + 1)th paddy field (WG(i) > d), the lateral water loss through the same bund from the ith paddy field, could go mostly to the (i + 1)th paddy field. In this case, the Dupuit equation can be applied to estimate Lout,through . Extending the hardpan level in the ith paddy field to the (i + 1)th paddy field (Fig. 5b), the lateral flow rate per unit length of the bund between the ith and (i − 1)th paddy fields, qthrough(i) , is given by (for WP > 0 and for i < n): qthrough(i) =
k 2 − [HPW(i+1) − Z(fz(i) − fz(i+1) )] } {H2 2w PW(i)
Lin,through(i) =
where k is the horizontal hydraulic conductivity of the bund, w the bund thickness, and HPW is the height of the water
AP(i)
Lout,through(i−1)
(11)
Rainfall received in a catchment above the toposequence goes mostly to deep percolation, while surface runoff goes to the paddy fields, so Lin,through is negligible at the top position in the toposequence (i.e., Lin,through = 0 for i = 1).
2.2.2.3. Lateral water flow over bunds between paddy fields (Lout,over and Lin,over ). Surface runoff (Lout,over ) occurs only when the potential depth of standing water on a given day exceeds the bund height. Assuming that the outflow process of surface runoff occurs after the processes of all other components of the water balance are completed, Lout,over is determined as follows: Lout,over = WP∗ + Winput − Woutput − h
(12)
where WP∗ is the free water level on the previous day, Winput the input water depth to the field (sum of R, Lin,through , Lin,over and Lrunon ), and Woutput is output water depth from the field (sum of ETpaddy , Dtotal , Lout,through ). When Lout,over value computed using Eq. (12) is negative, there is no surface runoff. Lout,over is regarded as input water to the paddy field below, so Lin,over is estimated by weighting Lout,over by the area ratio of paddy fields (for i > 1): Lin,over(i) =
2.2.3. (9b)
AP(i−1)
AP(i−1) AP(i)
Lout,over(i−1)
(13)
Total downward water flow (Dtotal )
The total downward water flow (Dtotal ) consists of downward water flow within the field (Dwithin ) and downward water flow under the bund (Dunder ). Both sub-components of Dtotal are empirically determined, as follows.
510
e c o l o g i c a l m o d e l l i n g 2 0 4 ( 2 0 0 7 ) 503–515
standing water is present above the soil surface (i.e., ε = 1 for WP ≥ 0).
Fig. 6 – Relationships between downward water flow and clay content of the topsoil when there is only one water table (the groundwater head is greater than or equal to the hardpan level) and when perched water is developwed in paddy fields (e.g., the groundwater table level of −500 mm).
2.2.3.1. Downward water flow within the field (Dwithin ). The estimation method for downward water flow through the hardpan within the field (Dwithin ) is based on data analysis methodology of field experiments in Northeast Thailand (Tsubo et al., 2005). Dwithin (mm day−1 ) is empirically related to clay content of the topsoil when no perched water is present in the field; Dwithin decreases as clay content of the topsoil increases (Tsubo et al., 2007). When perched water is present, Dwithin further increases as WG (mm) deepens, as follows: Dwithin =
8.05 − 0.03WG + 0.24 C
(14)
where C is clay content (%). This empirical equation was developed using the data collected from the experiments in Northeast Thailand. WG is likely to affect Dwithin until WG reaches a critical groundwater level (critical WG = −500 mm for the present study) (Tsubo et al., 2007). Beyond the critical WG , the groundwater table will not further affect Dwithin , so WG is assumed to be −500 mm for WG < −500 mm. The relationships between Dwithin and C for WG ≥ d and for WG = −500 mm (from Eq. (14)), are shown in Fig. 6. Downward water loss declines markedly when the standing water disappears and becomes zero as it reaches the hardpan level. Dwithin without the standing water is approximated using a parameter for reducing Dwithin , ε (0 ≤ ε < 1), which is determined from WP . When it is assumed that Dwithin is proportional to the WP for WP < 0, ε is defined as: ε=
d − WP d
(15)
The depth of standing water in a field does not have much affect on downward water movement, as standing water levels are relatively shallow throughout the growing season in the study region. Therefore, Dwithin remains constant when the
2.2.3.2. Downward water flow under the bund (Dunder ). Downward water flow under the bund (Dunder ) (the under-bund percolation loss (Tuong et al., 1994)), is significant with a deep groundwater table, when compared with lateral seepage loss (Walker and Rushton, 1984, 1986). An earlier study (Tsubo et al., 2005) showed that, on clayey soils in the lowland environment, there is generally only one water table throughout the growing season, implying that Dunder is very small. On sandy soils, a shallow perched water table may develop, especially in the higher toposequence positions (Tsubo et al., 2006). When the perched water develops (WG < d), the horizontal water flow into the two bunds parallel to the x-axis (2X/n), and the bund adjacent to the paddy field above (Yfy,centre ), is considered as Dunder . Dunder is estimated as follows: Dunder =
(2X/n) + Yfy,centre qunder AP
(16)
In the present study, empirical values for qunder were used. Tuong et al. (1994) for example, used 0.1 m3 day−1 for a clayey soil, while Walker and Rushton (1984) estimated it to be 0.144 m3 day−1 for a loamy soil. For sandy soils, qunder can be assumed to be greater than 0.17 m3 day−1 , based on the hydraulic conductivity mentioned above (k > 1 m day−1 ). In addition, when the standing water disappears, Dunder can be negligible (i.e., Dunder = 0 for WP ≤ 0).
2.2.4.
Paddy evapotranspiration (ETpaddy )
Evapotranspiration (ET), in general, depends on field water conditions and crop growth stage. The FAO crop coefficient (Kc ) approach (Allen et al., 1998) is generally considered the best among crop ET models (Bouman et al., 2001) and is used in the present model. Because wind speed data are not readily available in the study region, reference evapotranspiration (ETo ) was estimated using the Priestley and Taylor (1972) model. The Priestley and Taylor (1972) model is limited to humid areas and could therefore be valid for estimating ETo in the lowlands during the wet-season. Where solar radiation data are not available, solar radiation can be estimated from sunshine ˚ ¨ equation. Alternatively, ETo can be hours using the Angstr om estimated from pan evaporation (Allen et al., 1998). Evapotranspiration, in general, declines with decreasing plant available water in the soil. ETpaddy under water stressed conditions can be adjusted simply by a water stress factor. In non-rice field crops, the stress factor depends on the water content of the soil, and plants become stressed when the water content reaches a critical level of plant available water (Ritchie et al., 1972). In the lowlands, water evaporates from the free water surface when standing water is present in the field, while water evaporates from the soil surface when standing water disappears from the field. Rice plants do not appear to become stressed immediately after the standing water disappears; however, canopy transpiration rate is reduced with any decline in plant available water (Fukai, 1995). This implies that WP can be used to determine a coefficient for reducing ETpaddy when WP is below the soil surface. It is assumed that ETpaddy declines linearly with the
511
e c o l o g i c a l m o d e l l i n g 2 0 4 ( 2 0 0 7 ) 503–515
(50% of ) is used to approximate WP for WP < 0. The porosity can vary from 0.4 m3 m−3 for a sandy soil (Sharma et al., 1995; Tsubo et al., 2006) and 0.5 m3 m−3 for a clayey soil (Aimrun et al., 2004). For a loamy soil, may range between 0.4 and 0.5 m3 m−3 ; 0.45 m3 m−3 was used for the present study.
2.3.
Fig. 7 – Graphical expression of the relationship of the plant water availability coefficient (Ks ) with free water level (WP ) when standing water disappears from paddy fields (Ks,d : the water stress coefficient at the hardpan level).
drop of WP from the soil surface to the level of the hardpan. Therefore, the stress factor, Ks , is given by (Fig. 7): Ks = 1 −
1 − K
s,d
d
Water levels (WG and WP )
Groundwater table in a hillslope can be estimated by linear interpolation between the water levels at the top and bottom of the slope (Wang and Anderson, 1982). Based on this approximate approach, when the perched water develops in the toposequence, the groundwater table with respect to the soil surface for the ith position (WG(i) ), is given by:
WG(i) = WG(1) +
m
(17)
where Ks,d is Ks when free water level is at the hardpan level. When free water disappears from a top layer of the soil (the top 0.5 m for the present study), ETpaddy can be negligible. Therefore, Ks,d is defined as [1−(d/−0.5)] for the present study. The product of Ks and Kc ETo gives ETpaddy for WP < 0. When standing water is present in the field, ETpaddy is at its potential level: i.e., Ks = 1 for WP ≥ 0.
2.2.5.
Model inputs used for the test are listed in Table 1. The soil in the Kampong Cham toposequence was loamy, so k, qunder and used for the model validation were, 0.8 m day−1 0.14 m3 day−1 and 0.45 m3 m−3 , respectively, as indicated earlier. Each component of the water balance model was calculated on a daily basis from early August to late October 2004. WP calculated using the water balance equation was compared with WP measured at the 1st (top), 6th (middle-1), 7th (middle-2), 8th (middle-3) and 17th (bottom) positions of the toposequence (Fig. 2c). The agreement of measured (observed) and predicted (calculated) values of WP was statistically evaluated using the d-index (Willmott, 1981):
d=1− WP
WG(n) − WG(1) n−1
(i − 1)
(18)
where WG(1) and WG(n) , are groundwater levels at the top and bottom positions, respectively. In the present study, measured WG(1) and WG(n) , were used for the interpolation. Estimates of WP are separated mainly into two stages: it is either above or below the soil surface (WP ≥ 0 or WP < 0, respectively). WP can vary between h (the maximum water level) and d (the minimum water level). WP on a given day is calculated as WP on the previous day (WP∗ ) plus S from the previous day. The steps for computation of WP are as follows: (i) defining if WP∗ is above or below the soil surface; (ii) defining whether WP∗ + S is above or below the soil surface; (iii) whether WP∗ + S is beyond the upper or lower limits (h or d, respectively). When the standing water disappears, the soil water content above the free water surface is below paddy soil porosity (). For the case where WP is positioned below the soil surface (i.e., WP∗ + S < 0), mean soil water content
Model evaluation
j=1
m
j=1
(pi − oi )
2
(|pi | + |oi |)
2
(19)
where oi and pi are the observed and predicted data, m is the number of the paired set data, oi = oi − o¯ and pi = pi − o¯ , and o¯ is the measured mean. When d = 1, the predicted data perfectly matches the observed data. The paired data set for the period from 7 August to 23 October was used for the model test. The sensitivity analysis was carried out by varying one of the four parameters affecting field water loss: the hydraulic conductivity of the bund (k), the under-bund percolation rate (qunder ), the runoff coefficient (˛), and minimum rainfall required for runoff to occur (ˇ/˛). k and qunder values for clayey and sandy soils were chosen for the analysis. As mentioned above, k and qunder used for the sensitivity analysis were 0.5 m day−1 and 0.10 m3 day−1 , respectively, for clayey soils, and 1.0 m day−1 and 0.17 m3 day−1 for sandy soils. ˛ and ˇ/˛ for the sensitivity analysis were ±50% of ˛ (= 0.3) and ˇ/˛ (=10 mm) values used for the model validation.
3.
Results and discussion
At the bottom position in the toposequence, the water level measured in the groundwater tube was equivalent to that of surface water, indicating that there was only one water table throughout the season in the lower part of the toposequence (Figs. 8 and 9). For the top position, on the other hand, during the first half of the study period, early August to late October, the groundwater level was below the level of the hardpan. However, for the second half of the study period, there was generally only one water table, as heavy rain was received in the period from 16 September to 6 October (334 mm; 55% of the total rainfall for the experimental period). The standing water in the paddy fields disappeared temporarily in late August and again in mid September, as there was little rainfall during those periods; there was no standing water remaining about 2–3 weeks after rainfall ceased in early October.
512
e c o l o g i c a l m o d e l l i n g 2 0 4 ( 2 0 0 7 ) 503–515
Fig. 8 – Seasonal changes in measured and simulated groundwater levels at the five toposequence positions (top, middle 1, middle 2, middle 3 and bottom) in Kampong Cham, Cambodia in the 2004 season.
Fig. 9 – Seasonal changes in measured and simulated paddy field water levels at the five toposequence positions (top, middle 1, middle 2, middle 3 and bottom) in Kampong Cham, Cambodia in the 2004 season.
Groundwater levels were estimated for middle-1, middle-2 and middle-3 paddies from values measured at top and bottom positions. Estimated WG for middle-1 and middle-2 were in general agreement with measured WG during most of the 3month period, except for mid September (Fig. 8). However, WG estimated for middle-3 was lower throughout the period. For this toposequence position, Dwithin could be overestimated, but there was a small variation in Dwithin across the toposequence positions: for the 3-month period, Dwithin = 43 mm at bottom to 58 mm at top (Fig. 10). Also, the underestimation of ETpaddy might have been small on account of standing water being present in the field during most of the period (Fig. 11; ETpaddy = 410 mm for top and 422 mm for bottom positions). The errors associated with the WG estimates were minimal. However, a large error might have occurred in Dunder for the middle-3 position, on account of a large variation in Dunder across the toposequence positions (Fig. 10; Dunder = 222 mm for top and 17 mm for bottom). Measured WG for middle-3 was around the soil surface mostly during the period, but the simulated WG was much lower. So, simulated Dunder could be
largely overestimated because qunder was very dependent on the depth of WG . Although estimates of WG at the top and bottom toposequence positions were beyond the aim of the present study, the change in WG could be estimated from rainfall data (Pannangpetch, 1993). The variation of WP from mid-August and late September was well predicted for the top, middle-1 and middle-2 positions (Fig. 9). For the middle-3 position, the measured WP fluctuated greatly during the whole period, but the trend was similar to that estimated by the model. For the bottom position, a negative WP was measured in late August and mid September, but the model estimated a positive WP until the standing water disappeared permanently in later October. It might therefore be better to use greater qthrough(n=17) in simulating more water drainage from the bottom paddy field. Information on discharge rate from the toposequrnce can be helpful for determining qthrough at the bottom field. For the top and middle-3 positions, the estimated dates for the disappearance of standing water from the field agreed with the actual measurements, while for middle-1, middle-2 and bottom posi-
513
e c o l o g i c a l m o d e l l i n g 2 0 4 ( 2 0 0 7 ) 503–515
Fig. 10 – Relationships of seasonal downward water losses within the field and under the bund with paddy field positions of the toposequence in Kampong Cham, Cambodia in the 2004 season.
Fig. 11 – Seasonal water budget for the toposequence in Kampong Cham, Cambodia in the 2004 season. Positive and negative values of the net lateral water flow indicate water loss and gain, respectively.
tions, the water disappearance dates estimated were about 1 week later than the actual dates measured (Fig. 9). The d-index was greater than 0.6, as shown in Table 2. This indicates that the model performed reasonably well for all the toposequence positions, especially in the upper part of the toposequence. The model was very sensitive to qunder for the middle and bottom toposequence positions, and k at bottom, but neither ˛ nor ˇ/˛ affected the model output (Table 2). The model started computing the water balance from the top to bottom positions. As a result, the model simulation errors might increase from higher to lower toposequence positions because the lateral water flow (seepage through the bund and runoff over the bund) from a field becomes input to one below the field. In fact, the model was very sensitive to those two coefficients (qunder and k), which were used to calculate Dunder and Lout,through (largely variable components of the water balance).
The ratio of ETpaddy to R (ETpaddy /R) was 0.67 for top and 0.69 for bottom (Fig. 11). These ratios were similar to ETpaddy /R observed in Northeast Thailand (Tsubo et al., 2005) and higher than ETpaddy /R in Southern Laos (Tsubo et al., 2006). No toposequential effect was found on ETpaddy . The dominant component of the downward water movement was the under-bund percolation loss (Fig. 10). Dunder was large in the upper part of the toposequence, as perched water was often present during the season. As a result of the toposequential effect, Dtotal increased with increased elevation. In contrast, the seasonal Lnet decreased with increased elevation at the lower positions (Fig. 11). In the upper part of the toposequence with the exception of the top position, the paddy fields gained water through Lnet . The lateral water loss (sum of Lout,through and Lout,over ) decreased with increased elevation; the decline was particularly marked for the lower positions (Fig. 12).
Table 2 – The d-index [an index of the agreement of observed and calculated values (Willmott, 1981)] for the model validation and sensitivity analysis Top Model validation Sensitivity analysis Hydraulic conductivity of bund, k 0.5 m day−1 (for a clayey soil) 1.0 m day−1 (for a sandy soil)
Middle 1
Middle-2
Middle-3
Bottom
0.89
0.79
0.67
0.62
0.69
0.85 0.92
0.80 0.79
0.68 0.68
0.62 0.63
0.76 0.66
Under-bund percolation rate, qunder 0.10 m3 day−1 (for a clayey soil) 0.17 m3 day−1 (for a sandy soil)
0.85 0.92
0.77 0.82
0.60 0.72
0.60 0.64
0.61 0.77
Runoff coefficient (slope), ˛ 0.45 (−50%) 0.15 (+50%)
0.92 0.87
0.77 0.80
0.68 0.67
0.60 0.63
0.71 0.68
Minimum rainfall for runoff, ˇ/˛ 5 mm (−50%) 15 mm (+50%)
0.86 0.91
0.80 0.79
0.66 0.67
0.61 0.63
0.67 0.71
514
e c o l o g i c a l m o d e l l i n g 2 0 4 ( 2 0 0 7 ) 503–515
4.
Fig. 12 – Relationships of seasonal lateral water loss and gain with paddy field positions of the toposequence in Kampong Cham, Cambodia in the 2004 season.
A similar trend was found in the lateral water gain (sum of Lin,through , Lin,over and Lrunon ). These results indicate that more lateral water inflow and outflow occur in the lower part of the toposequence. Further, a large effect of the run-on from the catchment was found on the lateral water inflow and outflow process in the lower part of the toposequence (Fig. 12). The model indicated that: (i) there was not much difference in the net lateral seepage between the toposequence positions (except for top and bottom positions), but large amounts of water move from field to field through lateral seepage; (ii) surface runoff took place only for a few days during the study period in the lower part of the toposequence; however, the net amount of the runoff for each runoff event was large; (iii) the run-on from the catchment was larger in lower than higher positions in the toposequence (with the exception of the top position). These different toposequential effects generated the overall variation of Lnet across the toposequence positions. Comparing with existing detailed models for water movement in paddy fields, the model developed in the present study is less data demanding. For example, ORYZA2000 and FEMWATER require a coefficient related to hydraulic conductivity for each soil layer, grid or cell. The present model requires only two inputs k and Dunder throughout the simulation run. Those two values could be available from the previous studies, such as those cited above. Most of the field level models (e.g., ORYZA2000, CERES-Rice) are one-dimensional and cannot be employed for the case where toposequential effects on the water balance are large. Other existing models estimating water balance for large scale ricelands (e.g., a catchment area of 43.6 km2 in Becu et al., 2003) may not lend themselves easily to estimate the discrepancies in field water availability between upper and lower positions of a short-distance toposequence (less than 1 km). Thus, the model developed in the present study fulfills the need of assessing water availability for crop modelling in the relatively short-distance toposequence, which often exists in rainfed rice systems in the Mekong region.
Conclusions
A water balance model for rainfed lowland rice paddies on gently sloping land was developed. The model computed daily water gain and loss processes along a toposequence at a site in Cambodia, and performed reasonably well for estimates of water balance along the toposequence. Lateral water movement was a major factor governing the water balance in the lower part of the toposequence. In contrast, the dominant component for the upper part of the toposequence was downward water flow, particularly under-bund percolation loss as the perched water developed in the upper part of the toposequence. Sensitivity analysis showed that the model had some degree of sensitivity to the under-bund percolation rate (qunder ) and the hydraulic conductivity for the bund (k). The determination of these parameters could thus be crucial for the water balance modelling. Also, a careful determination of the lateral water flow rate (qthrough ) needs to be made in simulating water availability in the bottom paddy field. Furthermore, many assumptions which may be unique to the study area (Mekong region) were used in the development of this model. This model is therefore limited to the Mekong region.
Acknowledgements Financial support from the Australian Centre for International Agricultural Research, the University of Queensland, and the International Rice Research Institute, is greatly acknowledged. The authors are also indebted to Cambodian research collaborators who assisted with data collection in the project.
references
Aimrun, W., Amin, M.S.M., Eltaib, S.M., 2004. Effective porosity of paddy soils as an estimation of its saturated hydraulic conductivity. Geoderma 1221, 197–203. Allen, R.G., Pereira, L.S., Raes, D., Smith, M., 1998. Crop Evapotranspiration. Guidelines for Computing Crop Water Requirements. FAO Irrigation and Drainage Paper 56. FAO, Rome. Becu, N., Perez, P., Walker, A., Barreteau, O., Le Page, C., 2003. Agent based simulation of a small catchment water management in northern Thailand. Description of the CATCHSCAPE model. Ecol. Model. 170, 319–331. Bolton, F.R., Zandstra, H.G., 1981. A soil moisture-based yield model of wetland rainfed rice. IRPS No. 62. International Rice ˜ Research Institute, Los Banos, Philippines. Bouman, B.A.M., Kropp, M.J., Tuong, T.P., Wopereis, M.C.S., Ten Berge, H.F.M., van Laar, H.H., 2001. ORYZA 2000: Modeling Lowland Rice. International Rice Research Institute, Los ˜ Banos, Philippines. Chen, S.K., Liu, C.W., Huang, H.C., 2002. Analysis of water movement in paddy rice fields (II) simulation studies. J. Hydrol. 268, 259–271. Fukai, S., 1995. Crop physiological approaches to understanding rice production under water liming conditions. In: Ishii, R., Horie, T. (eds.), Crop Research in Asia: Achievements and Perspective, Proceedings of the 2nd Asian Crop Science Conference, 22–23 August 1995, Fukui, Japan, pp. 240–245. Fukai, S., Rajatsasereekul, S., Boonjung, H., Skulkhu, E., 1995. Simulation modeling to quantify the effect of drought for
e c o l o g i c a l m o d e l l i n g 2 0 4 ( 2 0 0 7 ) 503–515
rainfed lowland rice in Northeast Thailand. In: Proceedings of the International Rice Research Conference on Fragile Lives in ˜ Fragile Ecosystems, Los Banos, Laguna, Philippines, February 13–17, pp. 657–672. Fukai, S., Basnayake, J., Cooper, M., 2000. Modeling water availability, crop growth, and yield of rainfed lowland rice genotypes in northeast Thailand. In: Proceedings of the International Workshop on Characterizing and Understanding Rainfed Environments, Bali, Indonesia, December 5–9, 1999, pp. 111–130. IRRI, 2005. World Rice Statistics. International Rice Research Institute, Manila, Philippines.
. Mackill, D.J., Coffman, W.r., Garrity, D.P., 1996. Rainfed Lowland Rice Improvement. International Rice Research Institute, Manila, Philippines. McMennamy, J.A., O’Toole, J.C., 1983. RICEMOD: A physiologically based rice growth and yield model. IRRI Research Paper Series No. 87. IRRI, Manila, pp. 33. Mohmood, R., 1998. Air temperature variations and rice productivity in Bangladesh: a comprehensive study of the performance of the YIELD and the CERES-Rice models. Ecol. Model. 106, 201–212. Odhiambo, L.O., Murty, V.V.N., 1996. Modeling water balance components in relation to field layout in lowland paddy fields. I. Model development. Agric. Water Manage. 30, 185–199. Pannangpetch, K., 1993. Application of model simulation to evaluate rice production at the district level. In: Bouman, B.A.M., van Laar, H.H., Zhaoqian, W. (eds.), Agro-ecology of rice-based cropping systems, Proceedings of the International Workshop on Agro-Ecological Zonation of Rice. 14–17 April 1993, Hangzhou, P.R. China, pp. 16–26. Priestley, C.H.B., Taylor, R.J., 1972. On the assessment of surface heat flux and evaporation using large-scale parameters. Mon. Weather Rev. 100, 81–92. Reddi, L.N., 2003. Seepage in Soils: Principles and Applications. John Wiley and Sons, Hoboken, New Jersey. Ritchie, J.T., Burnett, E., Henderson, R.C., 1972. Dryland evaporative flux in a subhumid climate. III. Soil water influence. Agron. J. 64, 168–173.
515
Sharma, P.K., Ingram, K.T., Harnpichitvitaya, D., 1995. Subsoil compaction to improve water use efficiency and yields of rainfed lowland rice in coarse-textured soils. Soil Till. Res. 36, 33–44. Singh, U., Ritchie, J.T., Godwin, D.C., 1993. A User’s Guide to CERES-Rice, V2.10. International Fertilizer Development Center, P.O. Box 2040, Muscle Shoals, Alabama. Suzuki, K., Goto, A., Mizutani, M., Sriboonlue, V., 2003. Simulation model of rainfed rice production on sloping land in northeast Thailand. Paddy Water Environ. 1, 91–97. Tsubo, M., Fukai, S., Basnayake, J., Tuong, T.P., Bouman, B., Harnpichitvitaya, D., 2005. Estimating percolation and lateral water flow on sloping land in rainfed lowland rice ecosystem. Plant Prod. Sci. 8, 354–357. Tsubo, M., Basnayake, J., Fukai, S., Sihathep, V., Siyavong, P., Sipaseuth, Chanphengsay, M., 2006. Toposequential effects on water balance and productivity in rainfed lowland rice ecosystem in Southern Laos. Field Crops Res. 97, 209–220. Tsubo, M., Fukai, S., Basnayake, J., Tuong, T.P., Bouman, B., Harnpichitvitaya, D., 2007. Effects of soil clay content on water balance and productivity in rainfed lowland rice ecosystem in Northeast Thailand. Plant Prod. Sci. 10, 232–241. Tuong, T.P., Wopereis, M.C.S., Marquez, J.A., Kropff, M.J., 1994. Mechanisms and control of percolation losses in irrigated puddle rice fields. Soil Sci. Soc. Am. J. 58, 1794–1803. Wade, L.J., Fukai, S., Samson, B.K., Ali, A., Mazid, M.A., 1999. Rainfed lowland rice: physical environment and cultivar requirements. Field Crops Res. 64, 3–12. Walker, S.H., Rushton, K.R., 1984. Verification of lateral percolation losses from irrigated rice fields by a numerical model. J. Hydrol. 71, 335–351. Walker, S.H., Rushton, K.R., 1986. Water loses through the bunds of irrigated rice fields interpreted through an analogue model. Agric. Water Manage. 11, 59–73. Wang, H.F., Anderson, M.P., 1982. Introduction to Groundwater Modeling: Finite Difference and Finite Element Methods. Academic Press, San Diego. Willmott, C.J., 1981. On the validation of models. Phys. Geogr. 2, 184–194.