- Email: [email protected]
Evaluating sustainable groundwater management options using the MIKE SHE integrated hydrogeological modelling package
Environmental Modelling & Software 14 (1999) 129–140
Evaluating sustainable groundwater management options using the MIKE SHE integrated hydrogeological modelling package C. Demetriou *, J.F. Punthakey Centre for Natural Resources, Department of Land and Water Conservation (DLWC), PO Box 3720, Parramatta, NSW 2124, Australia Received 30 April 1998; accepted 31 July 1998
Abstract The Wakool Irrigation District consists of 74,000 ha of agricultural land, a large part of which is affected by a rising watertable. Irrigation in the district started in 1936 when the watertable depth was about 9 m. Over the past 30 years, the watertable has risen as a consequence of irrigation, wet winters and inadequate drainage. In addition, the Wakool Irrigation District is underlain by the Calivil and Renmark aquifers of the Murray Basin. Rising groundwater pressures in these deeper aquifers are also contributing to the rising shallow watertable. Since 1981, a sub-surface drainage scheme with 48 bores pumping shallow saline groundwater into a 2000 ha evaporation basin has been developed in two stages to control the watertable in the Wakool Irrigation District. This paper describes the development and application of an integrated hydrogeological model for the Wakool Irrigation District, where management of rising watertable levels and land salinisation is a continuing problem. The model development is considered to be an important part in the establishment of sustainable water management policies for the Wakool Irrigation District. The Wakool model was developed by using the MIKE SHE integrated catchment-modelling package. The developed model enables analysis of the complex hydrogeological regime in the region, and prediction of the environmental impacts of various management options. It is able to describe temporal and spatial variations in the exchange of water between the land surface, drainage and supply systems, and the aquifers within the area. Management options proposed in the Wakool Land and Water Management Plan have been analysed for the period between 1975 and 2020. Various scenarios such as the implementation of on-farm recycling ponds in conjunction with laser levelling, deep-rooted perennials, tree planting, installation of deep groundwater pumps and the effect of shallow groundwater pumping, were investigated. The results from these simulations indicate that the best option is the implementation of shallow pumping. 1999 Elsevier Science Ltd. All rights reserved. Keywords: Sustainable water management; Integrated hydrogeological modelling; MIKE SHE; Wakool; Soil salinisation
Software availability Name of software: Developer: Contact Address: Telephone:
MIKE SHE (Flow module ⫹ preprocessor) Danish Hydraulic Institute (DHI) Agern Alle 5, DK-2970 Horsholm, Denmark ⫹ 45 45 769555, Fax: ⫹ 45 45 762567,
Abbreviations: DLWC, Department of Land and Water Conservation; LAI, Leaf Area Index; LWMP, Land and Water Management Plan; PCG, Preconditioned Conjugate Gradient; SOR, Successive overrelaxation; WIDM, Wakool Irrigation District Model; WTSSD, Wakool Tullakool Sub-Surface Drainage * Corresponding author. Fax: ⫹ 61-2-9891-5884; e-mail: [email protected]
Year first available: Hardware required: Software required: Program language:
Program size: Availability and cost:
email: [email protected] 1996 (Windows version) IBM compatible PC Windows 95 or Win NT Main program (FORTRAN), front-end (C). Source code not provided. 80 Mb From DHI, $US18,500
1. Introduction A rising watertable, and subsequent land salinisation has been observed during the last three decades in many irrigation districts throughout the southwestern part of
1364-8152/99/$ - see front matter 1999 Elsevier Science Ltd. All rights reserved. PII: S 1 3 6 4 - 8 1 5 2 ( 9 8 ) 0 0 0 6 4 - 4
130
C. Demetriou, J.F. Punthakey / Environmental Modelling & Software 14 (1999) 129–140
New South Wales, Australia. Crop yield is significantly influenced by a shallow saline watertable lying within a critical depth of less than 2 m below the land surface. Due to the high salinity of the groundwater in the area, the required amount of irrigation water is transferred to the area from the Murray River via irrigation channels. Concern over increasing environmental problems has encouraged communities to develop a Land and Water Management Plan (LWMP) to ensure a sustainable agricultural industry. The various management options proposed in the LWMP depend on the individual characteristics of the districts, and include common features concerning surface drainage, infrastructure, floodplain management, sub-surface drainage schemes, and onfarm management. Several hydrological and hydrogeological studies have been carried out in order to increase our understanding of the hydrological systems and assess future conditions in the district (e.g. Bogoda et al., 1994). These include development of a comprehensive MIKE SHE based hydrological model (Storm and Punthakey, 1995), which proved to be a valuable tool not only in the planning phase, but also during implementation and operation of the selected options. The model described here is an enhancement of the model developed by Storm and Punthakey (1995) with more layers and better calibration. It covers an area of 3240 km2, has a spatial resolution of 2 ⫻ 2 km and is also based on the MIKE SHE integrated modelling system. This paper describes the development of the enhanced model and presents the results obtained from several scenarios. The model provides a complete description of the hydrogeological system in the Wakool Irrigation District. This system involves temporal and spatial variations in the exchange of water between the land surface, drainage and supply systems, and the aquifers within the area. The close interaction between evaporation and the shallow watertable is shown by the model as well as the effect of rainfall distribution and intensity on the rise of the watertable. Different scenarios on water management have been analysed for a period between 1975 and 2020. Results are presented on scenarios such as the implementation of on-farm recycling ponds in conjunction with laser levelling, deep-rooted perennials, tree planting, shallow groundwater pumping and deep groundwater pumping. Detailed description of the above watertable management options is given in Section 6.1.
2. Description of the area The Wakool Irrigation District is located in the central part of the Murray Darling River Basin at the southeastern part of NSW, Australia. The district is bounded by the Wakool River in the south and Edward River in the north (Fig. 1). Irrigation commenced in the Wakool Irri-
Fig. 1.
Location of the Wakool Irrigation District.
gation District in 1936 and expanded substantially to the current extent of approximately 220,000 ha. The climate in the district is semi-arid with an average rainfall of approximately 360 mm per year. Heaviest rainfalls occur mainly during the southern-hemisphere winter months of June to August, and a large part of the district is prone to occasional flooding. The topography is in general flat with a westerly direction of groundwater flow and natural surface drainage. Three major hydrogeological formations are located within the area. These are the Shepparton, Calivil and Renmark Formations (Bogoda et al., 1994). The Shepparton Formation constitutes the upper aquifer system. The depth and extent of the water bearing layers in the Shepparton are limited, with pronounced variations in transmissivity across the district. Prior streams and ancestral river systems, which are remnants of former riverbeds, play an important role in groundwater flow paths. The deep regional aquifer systems, the Calivil and Renmark Formations play an important role for the regional groundwater flow in the Murray Basin, and the piezometric heads in these aquifers control the groundwater exchange with the Shepparton Formation. Within the Wakool Irrigation District, there seems to be a fine balance between heads in the Shepparton and the Calivil/Renmark aquifers, which determines leakage between aquifers. In contrast to more easterly located irrigation areas, the piezometric heads in the deeper aquifers have reached the same overall level as the watertable, and therefore the potential for deep recharge within the area has decreased. In general, the Shepparton Formation provides a relatively poor drainage capacity for the infiltrating water through the rootzone. Thus, natural rainfall combined with large-scale irrigation in the district has caused the
C. Demetriou, J.F. Punthakey / Environmental Modelling & Software 14 (1999) 129–140
area of shallow watertable (less than 2 m from land surface) to increase from 7200 ha to 47,500 ha during the period of 1960 to 1975. A large-scale groundwater pumping scheme, the Wakool Tullakool Sub-Surface Drainage (WTSSD) system, which is shown in Fig. 1, was established in the early 1980’s and has partially alleviated the problems in the south-eastern part of the district. However, an area of 26,100 ha is currently affected by a high watertable, which is continuing to expand under the present irrigation management. The current irrigation practices and possible further rise in the piezometric head in the deeper aquifers are expected to cause a further spread of shallow watertable area in many parts of the district.
131
Fig. 2. Redistribution at the 2 km grid scale for the major crops within the model area.
3. Methodology The development of the Wakool model has involved the following steps: (1) identification of all relevant processes and information available; (2) data collation and processing; (3) setting up the model; (4) calibration of the model; and (5) performing model simulations to predict the effect of the proposed management options on the regional groundwater system. The hydrogeological information available for the Wakool Irrigation District is described in Table 1. For the simulation period, only two crop distributions, for 1987/88 and 1993/94, were available. The crop distribution information was obtained from satellite images. Due to the high cost associated with the purchase and analysis of satellite images, it was decided to utilise these two distributions only and assume that the 1987/88 was a representative distribution for the period between 1975 to 1992. The 1993/94 crop distribution (Fig. 2) was assumed to be representative of the crop distribution for simulation times from 1993 onwards. The crop distributions were then redistributed to the model’s 2 by 2 km grid. During the redistribution, the crop areas were con-
served, and also, crops were kept close to their original positions. The crops used in the model are (i) no crop (ii) rice, (iii) autumn crops and (iv) summer (perennial) crops. Lucerne and phalaris are examples of summer crops, while wheat and grazing oats are autumn crops. 3.1. Model conceptualisation A reliable groundwater mathematical model should reflect all the important natural processes and should provide the possibility to impose relevant human interventions on the system. Important features which are common for irrigation areas include: (1) natural and constructed drainage lines; (2) supply channels with particular emphasis on representation of sites where significant seepage losses may occur; and (3) groundwater pumping schemes, and pumping of saline groundwater to evaporation basins. The model must further be able to simulate the infiltration and evapotranspiration of the agricultural crops on specific soils as well as the capillary rise from shallow watertable areas. The groundwater system in the Wakool Irrigation Dis-
Table 1 Major information used in the model application Information (input)
Distribution
Comments
Rainfall and water usage
T and S
Potential evapotranspiration
T
Soil distribution Land use Hydrogeological characteristics Topography Supply channel network Drainage network Groundwater pumping Watertable observations
S S S S T and S S T and S T and S
Water usage for each crop type has been estimated for (1975–95) (pan evaporation) Obtained from the Deniliquin meteorological station situated 70 km away from Wakool Properties have been estimated from grain size distributions Obtained from remote sensing for year 1993/94 Well log information By interpolating the tops of well elevations The estimated seepage losses were imposed at selected locations within the model area Both natural and man-made channels within the farm area are included The pumping scheme commenced in 1981 Quarterly (sometimes half-yearly) data for approximately 1200 bores
T: temporal variations (daily estimates); S: spatial variations (averaged over the 2 km grid).
132
C. Demetriou, J.F. Punthakey / Environmental Modelling & Software 14 (1999) 129–140
trict area can be divided into two aquifer systems, Shepparton and Calivil/Renmark. The former is further divided into two layers, representing a sandy and relatively permeable layer in the upper 15 to 20 meters from the surface (Upper Shepparton Formation). It is assumed that the main lateral groundwater movement in the Shepparton Formation occurs in this layer. Beneath the sandy layer there is the low permeable Lower Shepparton Formation, which primarily is responsible for the vertical flow exchange between the Upper Shepparton and the Calivil/Renmark system. Within the Upper Shepparton Formation, the prior and ancestral streams form significant routes for groundwater flow. The available hydrogeological information on the Calivil/Renmark Formations is in general very sparse, and refinements in the calibration of the deeper aquifers will be required when additional information becomes available. However, since the flow exchange between these formations and the Shepparton Formation is an important factor, they have been included in the model. Thus, the developed model consists of three layers: (1) Upper Shepparton; (2) Lower Shepparton and; (3) Calivil/Renmark. 4. The numerical model The usefulness of the MIKE SHE (Danish Hydraulic Institute, 1993) was shown by previous work in New South Wales (Carr et al., 1993). It was chosen, for the Wakool model, over other packages such as the widely used MODFLOW (McDonald and Harbaugh, 1988) because MIKE SHE is able to model in more detail the interaction between evaporation and shallow watertable. It is also able to simulate in a more advanced way the interaction between surface water from channels and rivers and the groundwater system. MIKE SHE provides a much better evapotranspiration module with root extraction as well as a comprehensive module for open channel/river flow. In addition, it should be noted that MIKE SHE provides a coupled saturated and unsaturated zone while MODFLOW is not able to simulate the flow through the unsaturated zone. An unsuccessful attempt was made in the past by the authors to model an irrigation district with similar characteristics to those in the Wakool Irrigation District, utilising MODFLOW with its simplified version of evaporation. The calibration of the MODFLOW based model proved to be impossible and the model predictions were unreasonable. This experience indicated to us that due to the shallow watertable in the area, flow through the unsaturated zone is a very critical component in our model development, and it cannot be handled adequately by MODFLOW. 4.1. Components of MIKE SHE water movement module The Wakool model, described in this paper, is based on the MIKE SHE modelling package. MIKE SHE is an
integrated package and covers all the major processes of the hydrologic cycle. The MIKE SHE Water Movement module has been designed with a modular program structure comprising of different process-oriented components, each describing a major flow process of the hydrological cycle. Used in combination, they describe the entire hydrological cycle. The components utilised in the model are described in the following Sections 4.1.1–4.1.5. 4.1.1. Evapotranspiration and interception component The evapotranspiration component uses meteorological and vegetative input data to predict, with empirically based equations and on a spatially distributed basis, the total evapotranspiration and net rainfall amounts resulting from different processes which take place. These processes are interception of rainfall by the canopy, drainage from the canopy, evaporation from the canopy surface, evaporation from the soil surface, and finally, uptake of water by plant roots and its transpiration. This component interacts with the unsaturated zone component, providing net rainfall and evapotranspiration loss rates and using information on soil moisture conditions in the root zone. The model has been developed at the Royal Veterinary and Agricultural University in Denmark and published in Kristensen and Jensen (1975). The calculated actual evapotranspiration is based on potential evaporation rate and the actual soil moisture status in the root zone, which are required as input data. 4.1.2. Overland flow and channel flow component When the net rainfall rate exceeds the infiltration capacity of the soil, water is ponded on the ground surface. This water is available as surface runoff, to be routed down-gradient towards the river system. The exact route and quantity is determined by the topography and flow resistance as well as by losses due to evaporation and infiltration along the flow path. The water reaching the river system as surface and subsurface flow is routed downstream in a separate node-point system. Both the overland flow and the channel flow are modelled by approximations of the St. Venant equations (Ritzema, 1994). 4.1.3. Unsaturated zone component The unsaturated zone component links the two horizontal two- and three-dimensional surface and subsurface flow processes together. The flow is described by the one-dimensional governing equation, the so-called Richards’ equation (Ritzema, 1994). Soil water systems comprise three phases; solid, liquid and gas which introduce non-linear terms in the Richards’ equation. Knowledge about the soil physical properties is required in order to obtain a solution to Richards’ equation. The upper part of the unsaturated zone includes root extrac-
C. Demetriou, J.F. Punthakey / Environmental Modelling & Software 14 (1999) 129–140
tion for the transpiration process. This is explicitly incorporated in the equation by sink terms. The integral of the sinks over the entire root zone depth amounts to the total actual evapotranspiration. Soil evaporation is catered for in the first sink term below the land surface. The interaction between the unsaturated and saturated zone is solved by an iterative mass balance procedure. This coupling procedure ensures a realistic description of the watertable fluctuations in situations with shallow soils where accounting for a variable specific yield above the watertable is important. 4.1.4. Saturated zone component The saturated zone component of the MIKE SHE water movement module calculates the saturated subsurface flow in the catchment. MIKE SHE allows for fully three-dimensional flow in a heterogeneous aquifer with conditions shifting between unconfined and confined conditions. The spatial and temporal variations of the dependent variable (the hydraulic head) is described mathematically by the non-linear Boussinesq equation (Ritzema, 1994) and solved numerically by an iterative finite difference technique. MIKE SHE gives the opportunity to choose between two matrix solver modules: the SOR module based on a successive over-relaxation technique and the PCG groundwater module based on a preconditioned conjugate gradient solution technique. 4.1.5. River/aquifer exchange component A river system influences a large part of the groundwater system in a catchment as it traverses the catchment in many directions. The impact of the river on groundwater heads is in horizontal and vertical directions. The surface area of the river system is, however small compared to the catchment area, and for regional modelling the river width may typically occupy a small percentage of the grid cell size. In many applications, the river can therefore be represented in a separate node system running along the boundaries of the grid cells acting as a source/sink line. The exchange component assumes that the width of the river is small compared to the grid cell dimension. The river flow computations are carried out in the corners of the grid cells. Interaction between the river, the groundwater system and the overland flow is assumed to take place in the middle of the intermediate river links connecting adjacent computational nodes.
133
addition, the optimal timestep size varies in time as a consequence of different conditions in the hydrological regime within a catchment. During rather stationary (dry) periods the timesteps may be large compared to transient situations during wet periods. Provisions for allowing different timesteps in the different components have been made in MIKE SHE. The unsaturated zone, exchange and open channel overland flow components use identical timestep. MIKE SHE water movement module allows the user to specify larger timesteps for the saturated zone component since the results are less sensitive to timestep size in this component. 4.3. Application of MIKE SHE in the Wakool irrigation area MIKE SHE provides the scope to analyse a wide variety of management options. These options not only involve the groundwater system, but also changes on the land surface, such as: sealing of irrigation distribution system, extension of drainage network, on-farm practices etc. Due to unavailability of field data, flooding events from the Wakool and Edward rivers have not been introduced in the model, but the effects of local flooding have been analyzed. A finite difference grid with a spatial resolution of 4 km2 was used to represent the spatial variations of catchment characteristics, inputs and imposed stresses. The drainage system is described with the river module, whereas the supply system is represented only at reaches where seepage losses occur. In the Upper Shepparton Formation, groundwater flow direction is from the southeast to northwest (Fig. 3). Based on Fig. 3, it can be assumed that for modelling purposes, there is no flow across the Edward (north of the model) and Wakool (south of the model) rivers. Fixed head boundary conditions have been applied along the eastern and northwestern boundaries. For the rest of the model layers, fixed head boundaries were assumed
4.2. Timestep calculations The time scale for different flow processes in a catchment varies. For example, the surface flow component reacts much faster on a rainfall input than the subsurface flow component. Thus, optimal timestep size (i.e. largest possible timestep without introduction of numerical errors) is different for the individual components. In
Fig. 3. Boundary conditions and general contour levels indicating the direction of flow for the model area.
134
C. Demetriou, J.F. Punthakey / Environmental Modelling & Software 14 (1999) 129–140
for all the model boundaries. These assumptions are based on the water level bore records, which indicate that there is minimal change in heads across these boundaries. The model has been run with variable timesteps depending on the actual conditions. Because the model has to calculate overland flow and infiltration rates at the land surface, it was necessary to apply very short timesteps. The timestep size was adjusted automatically, and varied from a high value of 12 hours to as low as sixty seconds. To give an indication of computation time requirements, it took one hour of CPU time to simulate a six-year period on a Pentium 100 MHz computer.
5. Model calibration The model was calibrated over the period 1975 to 1995. This period includes both dry and wet years. No major flood impacts occurred, thereby providing a suitable period for the model calibration. Historical measurements of groundwater levels in more than 1200 bores screened in the sandy layers of the Shepparton Formation were used. Some data on streamflows and soil moisture was available for sufficient calibration of the channel and soil moisture components of the model. By including these processes in the model, we were able to find out if they play a major role in the water movement, and also to establish which data should be collected. The main emphasis of the calibration process has been to provide good temporal and spatial agreements between simulated and observed groundwater levels throughout the irrigation district. The calibration performance can be measured by the ability of the model to predict historical regional groundwater flow patterns and historic trends in groundwater levels over the calibration period. The results from the calibration runs are described in the following Sections 5.1–5.3.
Fig. 4. Comparison of observed and simulated groundwater levels at selected bores. The bore numbers are as they are stored in the database. They are located in Fig. 5 at coordinates (x,y,layer) as used in the model.
bore 287 over the period of 1975 to 1980 do not match the observed values. This can be attributed to flooding which, due to lack of information, was not taken into account in this model.
5.1. Graphs of groundwater levels vs time
5.2. Contour plots at selected times
The observed readings at different bores are compared with the ones computed by the model. Throughout the calibration process, a representative sample of 21 bores was selected for which the water levels were compared. Fig. 4 shows the time series of six such bores. The model results are shown as solid lines while the water level observations are shown as dots at the time the reading was taken. Fig. 4 shows that there is a general drop of the groundwater elevation which can be attributed to the groundwater pumping from the 48 pumps which have been operating since 1981. This observation is more obvious for the bores 230 and 389, which are situated within the area affected by the pumping scheme. From Fig. 4, it can be seen that the computed heads for
Fig. 5 presents a contour plot for May 1993. In Fig. 5 (a) the entire Wakool model is shown, while in Fig. 5 (b) a sub region of the model (Tullakool irrigation area) is presented in more detail. In Fig. 5, the evaporation basins are presented as well as the drainage channels. In order to evaluate the performance of the calibration of the model, plots similar to the one shown in Fig. 5 were plotted at approximately 5-year intervals. Since errors are accumulating over the calibration period, results at or near the end of the calibration period are expected to have maximum errors. For this reason, results at July 1994 (which is close to the end of the period of calibration) where selected to be presented in this paper. A calibration with an overall average error (difference
C. Demetriou, J.F. Punthakey / Environmental Modelling & Software 14 (1999) 129–140
135
Fig. 5. Comparison of observed (circles) and simulated groundwater heads at July 1994. (a) Model area, (b) Subsection of the model area showing in detail the evaporation basins (light blue), drainage channels (dark blue lines), extraction pumps (light blue dots). Colours show the groundwater levels as obtained from the model. The small circles indicated the observation bores, which were monitored at that particular time. A good calibration is achieved when there is a colour match between the colour in the circles and the colour of the surrounding area.
136
C. Demetriou, J.F. Punthakey / Environmental Modelling & Software 14 (1999) 129–140
between the computed results and observed values at any particular time at the observed points) of less than 1.0 m was achieved which can be considered as good taking into consideration the grid size of the model. Inspection of maps with simulated and observed water levels revealed that the model is able to achieve good agreement with the observed groundwater levels both spatially and temporally. It should be acknowledged that the developed model is based on various assumptions: boundary conditions; parameters for flow in the saturated and unsaturated zones; pumping stresses in the deeper aquifers; and others. The results for different scenarios produced by the model depend on these assumptions. Therefore, when more detailed information becomes available, the model should be recalibrated and the scenarios re-evaluated. 5.3. Percentage of the area affected by shallow watertable Another indicator, which was used in the model calibration, is the percentage of the area with watertable lying within 2.5 m from the ground surface. A comparison of the observed and computed percentage of the area affected by shallow watertable was carried out at July 1987 and 1995. The year 1987 is the year when the Wakool Sub-Surface Drainage (WSSD) scheme was fully operational and 1995 is the year when the Wakool LWMP commenced. Table 2 shows the computed results for these two landmark years as well as the corresponding observed values. As shown in Table 2, a high degree of agreement between the computed and the observed values has been achieved. 5.4. Sensitivity analysis of the Wakool model Sensitivity analysis was carried out to establish which parameters were most affecting the model. Model runs were carried out to investigate the model sensitivity on hydraulic conductivity, storage coefficient, pumping, leakage, and evapotranspiration.
Table 2 Comparison of the computed and observed percentage of the area affected by shallow watertable for July 1987 and 1995 % of area* in July 1987 Derivation of with watertable depths values (m BGL**)
Computed Observed *
% of area* in July 1995 with watertable depths (m BGL**)
⬍ 1.5
⬍ 2.0
⬍ 2.5
⬍ 1.5
⬍ 2.0
⬍ 2.5
9 9
13 12
18 16
5 6
7 7
12 13
The total model area is 810 (2 ⫻ 2 km) cells ⫽ 3240 km2; BGL: Below ground level.
**
5.4.1. Hydraulic conductivity The hydraulic conductivity was increased (from the spatially distributed calibrated values) by two orders of magnitude and decreased by two orders of magnitude but the change was in the order of a few centimeters up or down respectively which indicated that the effect was very small. This indicates that the process, which has the most pronounced effect on the watertable, is the recharge through the unsaturated zone. If the recharge through the unsaturated zone remains unchanged, the water level will not rise even with significant decrease in hydraulic conductivity. 5.4.2. Storage coefficient The storage coefficient for the confined layers was changed to values ranging from 0.5 to 1.25 times the calibrated spatially distributed values. The effect on the water level was insignificant (a few centimeters up or down) which indicates that the effect is minor. This again can be attributed to the fact that the recharge through the unsaturated zone is the process most affecting the watertable. 5.4.3. Pumping and channel leakage The pumping and channel leakage, were changed from 80% to 120% of their spatially distributed calibrated values. The effect on the water level was significant (up to 0.5 m change in general, and in some cases up to 1.0 m change) which indicates that the effect is locally significant. 5.4.4. Evapotranspiration The variation in evapotranspiration between 85% to 100% of the time-varying pan value had a significant impact on water levels (more than 1.0 m change). This indicates that this was a significant parameter. Field measurements of this parameter are needed in order to reduce the uncertainty in the calibrated model. Moreover, on-farm management options, which impact on evapotranspiration, may result in a significant change in the water levels. Based on the above descriptions, it can be said that the model is very sensitive to pumping, leakage and evapotranspiration. Pumping and leakage have a local effect while evapotranspiration has a regional impact.
6. Predictive simulations To determine the impact of changing management options over time, the timeframe for the model was extended by 25 years from 1995 to 2020. Before the commencement of the simulations, records (time series) on rainfall/irrigation, evapotranspiration, crop distribution, leaf area index, root depth and pumping/leakage need to be provided to the model up to the year 2020.
C. Demetriou, J.F. Punthakey / Environmental Modelling & Software 14 (1999) 129–140
These records are available up to 1995. Thus, artificial records were constructed to the year 2020. (a) Rainfall. The record was developed utilising the historical daily rainfall information at the Deniliquin meteorological station (70 km away from Wakool) extending from 1878 to 1996. A continuous 25-year portion of the record, which represented the average rainfall of the whole record, was established by utilising the moving average method with a moving average of 25 years. For an average rainfall over a 25-year period, the position of the average value in the moving average series indicates the relative position of the middle of the 25-year record in the actual rainfall record. This record was added on top of 1995 to extend the record to year 2020. A similar approach can be adopted for selecting a set of 25 “dry” years (by selecting the minimum value of the moving-average, or a set of 25 “wet” years (by selecting the maximum value of the moving average). (b) Irrigation. Based on the rainfall record and crop water demands, irrigation time series for each crop were developed. The crops used in the model are: no crop; rice; perennial (summer); and autumn crops. (c) Evapotranspiration. The evaporation data was obtained from the Deniliquin meteorological station, which has records from 1977 onwards. Thus, the evapotranspiration series for the selected 25-year period was synthesised using statistical methods based on the records of rainfall and evapotranspiration from 1978 to 1995. Knowing the evapotranspiration record, an actual evapotranspiration series for each crop was generated using the crop coefficients. (d) Crop distribution. The crop distribution of 1993/94 was used from 1993 onwards. The crop distribution was obtained from satellite images and consisted of the crops mentioned above. (e) Leaf area index and root depth. Both records were extended by repeating the existing yearly variations to year 2020. (f) Pumping/leakage records. In order to satisfy all the requirements by the different scenarios examined, various records for pumping and leakage rates were constructed. 6.1. Simulation of various management scenarios Various combinations of the following management options were simulated with the model. 6.1.1. Recycling ponds in conjunction with land forming (laser levelling) On-farm recycling ponds are specially constructed ponds in every farm, which will collect any surface runoff from rainfall and excess irrigation. For our modelling
137
exercise, the effect of the recycling ponds is simulated as recharge to the watertable caused by leakage. Due to the big grid size of the model, the effect of the recycling ponds is implemented in the model as points, which leak into the aquifer. The recharge rate was estimated to be 46 mm/year (Demetriou et al., 1998) and the area affected is proportional to the area of the farm. Laser levelling will be carried out in conjunction with the recycling ponds option. The aim of laser levelling is to contour the farm in such a way so that seepage will be minimised. Land forming (laser levelling) was implemented in the model by increasing the Manning surface number (Ritzema, 1994) in the model (the soil surface becomes smoother and therefore will be less seepage). The only information available on the implementation, was that the percentage of laser levelled area within the model, will start in 1995 and reach 45% at year 2000 and 75% at year 2005. Based on this, the Manning surface number was increased from 1 to 20 for randomly selected model grid cells in these two model periods. In addition, the detention storage coefficient (Ritzema, 1994) was modified. For the calibration period 1975 to 1995, it was set to 1.0 m for the cells covering the evaporation basins (see Fig. 2), 0.2 m for the rice cells and 0.01 m for the rest of the model area. From 1995 onwards, the detention storage coefficient for model cells other than those representing rice and evaporation basins, was decreased to 0.0075 m. 6.1.2. Deep rooted perennials The objective of planting deep-rooted perennials is to extract water through the root zone and release it to the atmosphere by transpiration. The effectiveness of this scheme is limited by the root-depth and the salt tolerance of the plant. For this model scenario, the root depth of the summer (perennial) crops was increased from 0.5 m to 2.0 m. This option becomes active in the model from 1995 onwards. 6.1.3. Tree planting The objective of tree planting is to capture, through the root zone, any freshwater which might leak from irrigation channels. The effectiveness of this scheme is depended on the density and depth of the roots as well as the salt tolerance of the trees. Line extractions were placed along the irrigation channels (100 m on either side) within the Tullakool Irrigation area operating at an extraction rate of 1.0 mm/day (200 m3/day/km—strip). The width of the strip is 200 m, 100 m on either side from the irrigation channels. The water extraction was dropped to 35% during the winter months. The trees are modeled as pumps, and they extract water over the whole depth of the top layer (Upper Shepparton Formation). Tree planting becomes active from 1995 onwards and the trees are assumed fully mature from the time of planting.
138
C. Demetriou, J.F. Punthakey / Environmental Modelling & Software 14 (1999) 129–140
6.1.4. Shallow pumping The objective of shallow groundwater pumping is to pump the salty groundwater to the evaporation basins and thus slow or reverse the rise of the watertable. An extra, 48 shallow pumps were added to the model and were installed within the Tullakool irrigation area. Their locations were nominated by the office of DLWC in Deniliquin. The pumping rate for each pump is 250 ⫻ 103 m3/year starting from 1995. 6.1.5. Deep groundwater pumping options The objective of deep groundwater pumping is to reverse the upward movement of salty water from the lower to the upper aquifer. The effect of deep groundwater pumping is expected to be a long term one. Three deep pumps (extract water from the Calivil/Renmark Formations) were installed within the Tullakool irrigation area at the locations nominated by the office of DLWC Deniliquin office. The rate is 7.3 ⫻ 106 m3/year and starts from 1995 onwards. 6.2. Results of scenario modelling Table 3 presents a summary of the results from the different scenarios examined. The percentage of the area with watertable lying within 1.5 m, 2.0 m and 2.5 m from the land surface is presented in Table 3. Scenario 1 shows the results for the “basic” (as in 1995) scenario for the years 2005 and 2020. At year 2020, the area with watertable within 2 m reaches 29% of the total area. Scenarios 2 to 6 present the results obtained from recycling pond in conjunction with other options. Under these scenarios, the best performance (the least percentage of the area with the watertable lying within
2.0 m of the ground surface) was achieved by scenario 6 (48-pump option). Scenarios 7 and 8 cover the cases where all the options are incorporated at different implementation rates and with different detention storage. In scenario 7 the detention storage coefficient reduced to 50% of the ‘as in 1995’ case, while in scenario 8 it was reduced to 75%. As expected, slightly better performance was achieved by the smaller detention storage. Scenarios 9 and 10 investigate the effect the rainfall variation has on the results. The results are very sensitive to the average annual rainfall but more importantly to rainfall intensity and distribution with respect to time. Scenarios 11 to 14 investigate the effect of the different options as if they were implemented individually. The best performance of these was obtained from the 48-pumps option (scenario 13). 7. Model audit by external reviewers External consultants reviewed the developed model and a report with comments and recommendations from the reviewers as well as from the modelling team, was produced (Meynink et al., 1997). The following are some of the recommendations from the review process: 쐌 Field measurements need to be carried out in order to establish better estimates for the important parameters. These parameters are: evapotranspiration and leaf area index; surface roughness coefficient; detention storage coefficient. For deep-rooted perennials, information is needed on root development, root density and tolerance to groundwater salinity, and for tree planting, capacity of the plant to extract water, rates of growth of tree and salt tolerance.
Table 3 Percentage of the area affected by shallow watertable under various simulated management scenarios for March 2005 and 2020 Simulated scenario
% of area* in March 2005 with watertable depths (m BGL**) ⬍ 1.5
1. Basic (as in 1995) 2. Rec. Ponds ⫹ (land forming) LF1 3. Rec. Ponds ⫹ (land forming) LF2 4. Rec. Ponds ⫹ LF1 ⫹ deep roots 5. Rec. Ponds ⫹ LF1 ⫹ trees 6. Rec. Ponds ⫹ LF1 ⫹ 48 pumps 7. Rec. Ponds ⫹ LF1 ⫹ 48 pumps ⫹ 3 deep pumps ⫹ deep roots 8. Rec. Ponds ⫹ LF2 ⫹ 48 pumps ⫹ 3 deep pumps ⫹ deep roots 9. Rec. Ponds ⫹ LF2 ⫹ 48 pumps ⫹ 3 deep pumps ⫹ deep roots (dry) 10. Rec. Ponds ⫹ LF2 ⫹ 48 pumps ⫹ 3 deep pumps ⫹ deep roots (wet) 11. Basic (as in 1995) ⫹ deep roots 12. Basic (as in 1995) ⫹ trees 13. Basic (as in 1995) ⫹ 48 pumps 14. Basic (as in 1995) ⫹ 3 deep pumps
5 3 4 3 2 2 2 2 1 3 5 3 3 3
⬍ 2.0 9 6 7 6 5 3 2 2 2 7 7 5 4 5
⬍ 2.5 12 10 11 10 9 7 5 5 4 16 12 10 9 10
% of area* in March 2020 with watertable depths (m BGL**) ⬍ 1.5 20 7 10 7 4 2 2 2 1 2 13 9 1 10
⬍ 2.0 29 10 14 10 8 4 3 3 1 4 23 17 17 19
⬍ 2.5 40 16 22 14 12 8 3 8 2 10 34 23 27 28
The total model area is 810 (2 ⫻ 2 km) cells ⫽ 3240 km2; LF1: Land forming with detention storage coefficient reduced by 50% of the "as in 1995" case; LF2: Land forming with detention storage coefficient reduced by 75% of the "as in 1995" case; **BGL: Below ground level.
*
C. Demetriou, J.F. Punthakey / Environmental Modelling & Software 14 (1999) 129–140
쐌 Continue water level and water quality monitoring, preferably at quarterly intervals, or, at least take readings twice a year. 쐌 Expand the deep bore network and continue water level and water quality monitoring (minimum twice a year, preferably quarterly). 쐌 Field verifications of the leakage locations and rates need to be conducted. 쐌 In order to obtain an indication on the validity of model outputs, uncertainty analysis should be carried out. This requires a great amount of extra work because many different models need to be set up and calibrated.
139
From the results obtained, it can be concluded that the model is very sensitive to parameters such as evapotranspiration, rainfall pattern, and the parameters associated with land forming (Manning surface number and the detention storage). The results are strongly influenced by the above parameters and field measurements need to be carried out in order to establish better estimates for them. In addition, the rise and fall of the watertable is strongly influenced by the intensity, distribution and duration of the rainfall. Despite the data limitations identified in Section 7, the use of MIKE SHE to develop a model for the Wakool Irrigation District has provided an important tool for evaluating possible future strategies for watertable reduction in this area.
8. Conclusions and recommendations The MIKE SHE modelling package was chosen over other packages such as the widely used MODFLOW, because it is well suited to describe the dynamic interaction between the surface and subsurface water systems, which is particularly important in the high watertable areas of the Wakool Irrigation District. A range of groundwater management options proposed in the land and water management plan (LWMP) for the Wakool Irrigation District have been investigated and their effectiveness have been evaluated. These options range from revegetation-type alternatives such as planting trees, salt tolerant deep rooted perennials, to engineering alternatives such as pumping from different aquifers, laser levelling (land forming), and utilising recycling ponds. The model predicts that with the current irrigation practices, the area with a watertable less than 2 m deep will further expand, and by the year 2020, will cover about 29% of the Wakool Irrigation District, an increase of 22% compared to that at 1995. Of all the groundwater management options assessed in this study, on-farm recycling in conjunction with land levelling (scenario 3) and the 48-pump option (scenario 13) perform the best. Under these two scenarios, the shallow watertable, at year 2020, is expected to reach 14% and 17% of the area respectively (an increase of 7% and 10% compared to 1995). It should be pointed out that the confidence in the 48-pump option (scenario 13) is very high because the pumping rates are known. On the other hand, the results from the on-farm recycling option (scenarios 2 and 3) carry a high degree of uncertainty since the coefficient for surface roughness and detention storage are not well known. Also, high uncertainty is associated with the deep-rooted perennials option (scenarios 4 and 12). This is due to many factors associated with the root development such as root density, root depth and plant tolerance to groundwater salinity.
Acknowledgements Mr. Borge Storm, consultant, Lawson and Treloar Pty Ltd, is acknowledged for his input in setting up the model and providing ideas on implementing and analysing the scenario simulations. Mr. K.R. Bogoda and Mr. N. Kulatunga, DLWC, Deniliquin are thanked for their advice and provision of field data. Mr. H. Schroo, DLWC, Parramatta is acknowledged for his constructive criticism and input during the model development and analysis of the scenario simulations.
References Bogoda, K.R., Kulatunga, N., Hehir, K., 1994. Overview of hydrogeology and assessment of subsurface drainage option for watertable control in Wakool. Internal Working Document, DLWC, pp. 123. Carr, R.S., Punthakey, J.F., Cooke, R., Storm, B., 1993. Large-scale catchment simulation using the MIKE-SHE model: 1. Process simulation of an irrigation district. International Conference on Environmental Management. Wollongong, Australia, Balkema, Rotterdam, pp. 467–472. Danish Hydraulic Institute, 1993. MIKE SHE water movement-user’s guide and technical manual, ed. 1.0. DHI, Denmark, pp. 81. Demetriou, C., Storm, B., Punthakey, J.F., 1998. Development of an integrated model to simulate on-farm water recycling in the Wakool Irrigation District. Second International Conference on Environmental Management (ICEM2), vol. 1. Wollongong, Australia, pp. 633–640. Kristensen, K.J., Jensen, S.E., 1975. A model for estimating actual evapotranspiration from potential evapotranspiration. Nordic Hydrology 6, 70–88. McDonald, M.G., Harbaugh, A.W., 1988. A modular three-dimensional finite-difference ground-water flow model. U.S. Geological Survey techniques of water resources investigations, Book 6, pp. 586. Meynink, W.J.C., Doherty, J., Demetriou, C., Punthakey, J.F., 1997. NSW DLWC Wakool SHE model review. Report # psm— wakool—051.1.doc, Pells Sullivan Meynink (Brisbane). Ritzema, H.P., 1994. Drainage principles and applications. International Institute for Land Reclamation and Improvement (ILRI) publication 16, 2nd ed. The Netherlands, pp. 1125.
140
C. Demetriou, J.F. Punthakey / Environmental Modelling & Software 14 (1999) 129–140
Storm, B., Punthakey, J.F., 1995. Modelling Environmental changes in the Wakool Irrigation District. MODSIM 95, Int. Congress on
Modelling and Simulation Proceedings, vol. 3. Newcastle, Australia, pp. 47–52.
Evaluating sustainable groundwater management options using the MIKE SHE integrated hydrogeological modelling package C. Demetriou *, J.F. Punthakey Centre for Natural Resources, Department of Land and Water Conservation (DLWC), PO Box 3720, Parramatta, NSW 2124, Australia Received 30 April 1998; accepted 31 July 1998
Abstract The Wakool Irrigation District consists of 74,000 ha of agricultural land, a large part of which is affected by a rising watertable. Irrigation in the district started in 1936 when the watertable depth was about 9 m. Over the past 30 years, the watertable has risen as a consequence of irrigation, wet winters and inadequate drainage. In addition, the Wakool Irrigation District is underlain by the Calivil and Renmark aquifers of the Murray Basin. Rising groundwater pressures in these deeper aquifers are also contributing to the rising shallow watertable. Since 1981, a sub-surface drainage scheme with 48 bores pumping shallow saline groundwater into a 2000 ha evaporation basin has been developed in two stages to control the watertable in the Wakool Irrigation District. This paper describes the development and application of an integrated hydrogeological model for the Wakool Irrigation District, where management of rising watertable levels and land salinisation is a continuing problem. The model development is considered to be an important part in the establishment of sustainable water management policies for the Wakool Irrigation District. The Wakool model was developed by using the MIKE SHE integrated catchment-modelling package. The developed model enables analysis of the complex hydrogeological regime in the region, and prediction of the environmental impacts of various management options. It is able to describe temporal and spatial variations in the exchange of water between the land surface, drainage and supply systems, and the aquifers within the area. Management options proposed in the Wakool Land and Water Management Plan have been analysed for the period between 1975 and 2020. Various scenarios such as the implementation of on-farm recycling ponds in conjunction with laser levelling, deep-rooted perennials, tree planting, installation of deep groundwater pumps and the effect of shallow groundwater pumping, were investigated. The results from these simulations indicate that the best option is the implementation of shallow pumping. 1999 Elsevier Science Ltd. All rights reserved. Keywords: Sustainable water management; Integrated hydrogeological modelling; MIKE SHE; Wakool; Soil salinisation
Software availability Name of software: Developer: Contact Address: Telephone:
MIKE SHE (Flow module ⫹ preprocessor) Danish Hydraulic Institute (DHI) Agern Alle 5, DK-2970 Horsholm, Denmark ⫹ 45 45 769555, Fax: ⫹ 45 45 762567,
Abbreviations: DLWC, Department of Land and Water Conservation; LAI, Leaf Area Index; LWMP, Land and Water Management Plan; PCG, Preconditioned Conjugate Gradient; SOR, Successive overrelaxation; WIDM, Wakool Irrigation District Model; WTSSD, Wakool Tullakool Sub-Surface Drainage * Corresponding author. Fax: ⫹ 61-2-9891-5884; e-mail: [email protected]
Year first available: Hardware required: Software required: Program language:
Program size: Availability and cost:
email: [email protected] 1996 (Windows version) IBM compatible PC Windows 95 or Win NT Main program (FORTRAN), front-end (C). Source code not provided. 80 Mb From DHI, $US18,500
1. Introduction A rising watertable, and subsequent land salinisation has been observed during the last three decades in many irrigation districts throughout the southwestern part of
1364-8152/99/$ - see front matter 1999 Elsevier Science Ltd. All rights reserved. PII: S 1 3 6 4 - 8 1 5 2 ( 9 8 ) 0 0 0 6 4 - 4
130
C. Demetriou, J.F. Punthakey / Environmental Modelling & Software 14 (1999) 129–140
New South Wales, Australia. Crop yield is significantly influenced by a shallow saline watertable lying within a critical depth of less than 2 m below the land surface. Due to the high salinity of the groundwater in the area, the required amount of irrigation water is transferred to the area from the Murray River via irrigation channels. Concern over increasing environmental problems has encouraged communities to develop a Land and Water Management Plan (LWMP) to ensure a sustainable agricultural industry. The various management options proposed in the LWMP depend on the individual characteristics of the districts, and include common features concerning surface drainage, infrastructure, floodplain management, sub-surface drainage schemes, and onfarm management. Several hydrological and hydrogeological studies have been carried out in order to increase our understanding of the hydrological systems and assess future conditions in the district (e.g. Bogoda et al., 1994). These include development of a comprehensive MIKE SHE based hydrological model (Storm and Punthakey, 1995), which proved to be a valuable tool not only in the planning phase, but also during implementation and operation of the selected options. The model described here is an enhancement of the model developed by Storm and Punthakey (1995) with more layers and better calibration. It covers an area of 3240 km2, has a spatial resolution of 2 ⫻ 2 km and is also based on the MIKE SHE integrated modelling system. This paper describes the development of the enhanced model and presents the results obtained from several scenarios. The model provides a complete description of the hydrogeological system in the Wakool Irrigation District. This system involves temporal and spatial variations in the exchange of water between the land surface, drainage and supply systems, and the aquifers within the area. The close interaction between evaporation and the shallow watertable is shown by the model as well as the effect of rainfall distribution and intensity on the rise of the watertable. Different scenarios on water management have been analysed for a period between 1975 and 2020. Results are presented on scenarios such as the implementation of on-farm recycling ponds in conjunction with laser levelling, deep-rooted perennials, tree planting, shallow groundwater pumping and deep groundwater pumping. Detailed description of the above watertable management options is given in Section 6.1.
2. Description of the area The Wakool Irrigation District is located in the central part of the Murray Darling River Basin at the southeastern part of NSW, Australia. The district is bounded by the Wakool River in the south and Edward River in the north (Fig. 1). Irrigation commenced in the Wakool Irri-
Fig. 1.
Location of the Wakool Irrigation District.
gation District in 1936 and expanded substantially to the current extent of approximately 220,000 ha. The climate in the district is semi-arid with an average rainfall of approximately 360 mm per year. Heaviest rainfalls occur mainly during the southern-hemisphere winter months of June to August, and a large part of the district is prone to occasional flooding. The topography is in general flat with a westerly direction of groundwater flow and natural surface drainage. Three major hydrogeological formations are located within the area. These are the Shepparton, Calivil and Renmark Formations (Bogoda et al., 1994). The Shepparton Formation constitutes the upper aquifer system. The depth and extent of the water bearing layers in the Shepparton are limited, with pronounced variations in transmissivity across the district. Prior streams and ancestral river systems, which are remnants of former riverbeds, play an important role in groundwater flow paths. The deep regional aquifer systems, the Calivil and Renmark Formations play an important role for the regional groundwater flow in the Murray Basin, and the piezometric heads in these aquifers control the groundwater exchange with the Shepparton Formation. Within the Wakool Irrigation District, there seems to be a fine balance between heads in the Shepparton and the Calivil/Renmark aquifers, which determines leakage between aquifers. In contrast to more easterly located irrigation areas, the piezometric heads in the deeper aquifers have reached the same overall level as the watertable, and therefore the potential for deep recharge within the area has decreased. In general, the Shepparton Formation provides a relatively poor drainage capacity for the infiltrating water through the rootzone. Thus, natural rainfall combined with large-scale irrigation in the district has caused the
C. Demetriou, J.F. Punthakey / Environmental Modelling & Software 14 (1999) 129–140
area of shallow watertable (less than 2 m from land surface) to increase from 7200 ha to 47,500 ha during the period of 1960 to 1975. A large-scale groundwater pumping scheme, the Wakool Tullakool Sub-Surface Drainage (WTSSD) system, which is shown in Fig. 1, was established in the early 1980’s and has partially alleviated the problems in the south-eastern part of the district. However, an area of 26,100 ha is currently affected by a high watertable, which is continuing to expand under the present irrigation management. The current irrigation practices and possible further rise in the piezometric head in the deeper aquifers are expected to cause a further spread of shallow watertable area in many parts of the district.
131
Fig. 2. Redistribution at the 2 km grid scale for the major crops within the model area.
3. Methodology The development of the Wakool model has involved the following steps: (1) identification of all relevant processes and information available; (2) data collation and processing; (3) setting up the model; (4) calibration of the model; and (5) performing model simulations to predict the effect of the proposed management options on the regional groundwater system. The hydrogeological information available for the Wakool Irrigation District is described in Table 1. For the simulation period, only two crop distributions, for 1987/88 and 1993/94, were available. The crop distribution information was obtained from satellite images. Due to the high cost associated with the purchase and analysis of satellite images, it was decided to utilise these two distributions only and assume that the 1987/88 was a representative distribution for the period between 1975 to 1992. The 1993/94 crop distribution (Fig. 2) was assumed to be representative of the crop distribution for simulation times from 1993 onwards. The crop distributions were then redistributed to the model’s 2 by 2 km grid. During the redistribution, the crop areas were con-
served, and also, crops were kept close to their original positions. The crops used in the model are (i) no crop (ii) rice, (iii) autumn crops and (iv) summer (perennial) crops. Lucerne and phalaris are examples of summer crops, while wheat and grazing oats are autumn crops. 3.1. Model conceptualisation A reliable groundwater mathematical model should reflect all the important natural processes and should provide the possibility to impose relevant human interventions on the system. Important features which are common for irrigation areas include: (1) natural and constructed drainage lines; (2) supply channels with particular emphasis on representation of sites where significant seepage losses may occur; and (3) groundwater pumping schemes, and pumping of saline groundwater to evaporation basins. The model must further be able to simulate the infiltration and evapotranspiration of the agricultural crops on specific soils as well as the capillary rise from shallow watertable areas. The groundwater system in the Wakool Irrigation Dis-
Table 1 Major information used in the model application Information (input)
Distribution
Comments
Rainfall and water usage
T and S
Potential evapotranspiration
T
Soil distribution Land use Hydrogeological characteristics Topography Supply channel network Drainage network Groundwater pumping Watertable observations
S S S S T and S S T and S T and S
Water usage for each crop type has been estimated for (1975–95) (pan evaporation) Obtained from the Deniliquin meteorological station situated 70 km away from Wakool Properties have been estimated from grain size distributions Obtained from remote sensing for year 1993/94 Well log information By interpolating the tops of well elevations The estimated seepage losses were imposed at selected locations within the model area Both natural and man-made channels within the farm area are included The pumping scheme commenced in 1981 Quarterly (sometimes half-yearly) data for approximately 1200 bores
T: temporal variations (daily estimates); S: spatial variations (averaged over the 2 km grid).
132
C. Demetriou, J.F. Punthakey / Environmental Modelling & Software 14 (1999) 129–140
trict area can be divided into two aquifer systems, Shepparton and Calivil/Renmark. The former is further divided into two layers, representing a sandy and relatively permeable layer in the upper 15 to 20 meters from the surface (Upper Shepparton Formation). It is assumed that the main lateral groundwater movement in the Shepparton Formation occurs in this layer. Beneath the sandy layer there is the low permeable Lower Shepparton Formation, which primarily is responsible for the vertical flow exchange between the Upper Shepparton and the Calivil/Renmark system. Within the Upper Shepparton Formation, the prior and ancestral streams form significant routes for groundwater flow. The available hydrogeological information on the Calivil/Renmark Formations is in general very sparse, and refinements in the calibration of the deeper aquifers will be required when additional information becomes available. However, since the flow exchange between these formations and the Shepparton Formation is an important factor, they have been included in the model. Thus, the developed model consists of three layers: (1) Upper Shepparton; (2) Lower Shepparton and; (3) Calivil/Renmark. 4. The numerical model The usefulness of the MIKE SHE (Danish Hydraulic Institute, 1993) was shown by previous work in New South Wales (Carr et al., 1993). It was chosen, for the Wakool model, over other packages such as the widely used MODFLOW (McDonald and Harbaugh, 1988) because MIKE SHE is able to model in more detail the interaction between evaporation and shallow watertable. It is also able to simulate in a more advanced way the interaction between surface water from channels and rivers and the groundwater system. MIKE SHE provides a much better evapotranspiration module with root extraction as well as a comprehensive module for open channel/river flow. In addition, it should be noted that MIKE SHE provides a coupled saturated and unsaturated zone while MODFLOW is not able to simulate the flow through the unsaturated zone. An unsuccessful attempt was made in the past by the authors to model an irrigation district with similar characteristics to those in the Wakool Irrigation District, utilising MODFLOW with its simplified version of evaporation. The calibration of the MODFLOW based model proved to be impossible and the model predictions were unreasonable. This experience indicated to us that due to the shallow watertable in the area, flow through the unsaturated zone is a very critical component in our model development, and it cannot be handled adequately by MODFLOW. 4.1. Components of MIKE SHE water movement module The Wakool model, described in this paper, is based on the MIKE SHE modelling package. MIKE SHE is an
integrated package and covers all the major processes of the hydrologic cycle. The MIKE SHE Water Movement module has been designed with a modular program structure comprising of different process-oriented components, each describing a major flow process of the hydrological cycle. Used in combination, they describe the entire hydrological cycle. The components utilised in the model are described in the following Sections 4.1.1–4.1.5. 4.1.1. Evapotranspiration and interception component The evapotranspiration component uses meteorological and vegetative input data to predict, with empirically based equations and on a spatially distributed basis, the total evapotranspiration and net rainfall amounts resulting from different processes which take place. These processes are interception of rainfall by the canopy, drainage from the canopy, evaporation from the canopy surface, evaporation from the soil surface, and finally, uptake of water by plant roots and its transpiration. This component interacts with the unsaturated zone component, providing net rainfall and evapotranspiration loss rates and using information on soil moisture conditions in the root zone. The model has been developed at the Royal Veterinary and Agricultural University in Denmark and published in Kristensen and Jensen (1975). The calculated actual evapotranspiration is based on potential evaporation rate and the actual soil moisture status in the root zone, which are required as input data. 4.1.2. Overland flow and channel flow component When the net rainfall rate exceeds the infiltration capacity of the soil, water is ponded on the ground surface. This water is available as surface runoff, to be routed down-gradient towards the river system. The exact route and quantity is determined by the topography and flow resistance as well as by losses due to evaporation and infiltration along the flow path. The water reaching the river system as surface and subsurface flow is routed downstream in a separate node-point system. Both the overland flow and the channel flow are modelled by approximations of the St. Venant equations (Ritzema, 1994). 4.1.3. Unsaturated zone component The unsaturated zone component links the two horizontal two- and three-dimensional surface and subsurface flow processes together. The flow is described by the one-dimensional governing equation, the so-called Richards’ equation (Ritzema, 1994). Soil water systems comprise three phases; solid, liquid and gas which introduce non-linear terms in the Richards’ equation. Knowledge about the soil physical properties is required in order to obtain a solution to Richards’ equation. The upper part of the unsaturated zone includes root extrac-
C. Demetriou, J.F. Punthakey / Environmental Modelling & Software 14 (1999) 129–140
tion for the transpiration process. This is explicitly incorporated in the equation by sink terms. The integral of the sinks over the entire root zone depth amounts to the total actual evapotranspiration. Soil evaporation is catered for in the first sink term below the land surface. The interaction between the unsaturated and saturated zone is solved by an iterative mass balance procedure. This coupling procedure ensures a realistic description of the watertable fluctuations in situations with shallow soils where accounting for a variable specific yield above the watertable is important. 4.1.4. Saturated zone component The saturated zone component of the MIKE SHE water movement module calculates the saturated subsurface flow in the catchment. MIKE SHE allows for fully three-dimensional flow in a heterogeneous aquifer with conditions shifting between unconfined and confined conditions. The spatial and temporal variations of the dependent variable (the hydraulic head) is described mathematically by the non-linear Boussinesq equation (Ritzema, 1994) and solved numerically by an iterative finite difference technique. MIKE SHE gives the opportunity to choose between two matrix solver modules: the SOR module based on a successive over-relaxation technique and the PCG groundwater module based on a preconditioned conjugate gradient solution technique. 4.1.5. River/aquifer exchange component A river system influences a large part of the groundwater system in a catchment as it traverses the catchment in many directions. The impact of the river on groundwater heads is in horizontal and vertical directions. The surface area of the river system is, however small compared to the catchment area, and for regional modelling the river width may typically occupy a small percentage of the grid cell size. In many applications, the river can therefore be represented in a separate node system running along the boundaries of the grid cells acting as a source/sink line. The exchange component assumes that the width of the river is small compared to the grid cell dimension. The river flow computations are carried out in the corners of the grid cells. Interaction between the river, the groundwater system and the overland flow is assumed to take place in the middle of the intermediate river links connecting adjacent computational nodes.
133
addition, the optimal timestep size varies in time as a consequence of different conditions in the hydrological regime within a catchment. During rather stationary (dry) periods the timesteps may be large compared to transient situations during wet periods. Provisions for allowing different timesteps in the different components have been made in MIKE SHE. The unsaturated zone, exchange and open channel overland flow components use identical timestep. MIKE SHE water movement module allows the user to specify larger timesteps for the saturated zone component since the results are less sensitive to timestep size in this component. 4.3. Application of MIKE SHE in the Wakool irrigation area MIKE SHE provides the scope to analyse a wide variety of management options. These options not only involve the groundwater system, but also changes on the land surface, such as: sealing of irrigation distribution system, extension of drainage network, on-farm practices etc. Due to unavailability of field data, flooding events from the Wakool and Edward rivers have not been introduced in the model, but the effects of local flooding have been analyzed. A finite difference grid with a spatial resolution of 4 km2 was used to represent the spatial variations of catchment characteristics, inputs and imposed stresses. The drainage system is described with the river module, whereas the supply system is represented only at reaches where seepage losses occur. In the Upper Shepparton Formation, groundwater flow direction is from the southeast to northwest (Fig. 3). Based on Fig. 3, it can be assumed that for modelling purposes, there is no flow across the Edward (north of the model) and Wakool (south of the model) rivers. Fixed head boundary conditions have been applied along the eastern and northwestern boundaries. For the rest of the model layers, fixed head boundaries were assumed
4.2. Timestep calculations The time scale for different flow processes in a catchment varies. For example, the surface flow component reacts much faster on a rainfall input than the subsurface flow component. Thus, optimal timestep size (i.e. largest possible timestep without introduction of numerical errors) is different for the individual components. In
Fig. 3. Boundary conditions and general contour levels indicating the direction of flow for the model area.
134
C. Demetriou, J.F. Punthakey / Environmental Modelling & Software 14 (1999) 129–140
for all the model boundaries. These assumptions are based on the water level bore records, which indicate that there is minimal change in heads across these boundaries. The model has been run with variable timesteps depending on the actual conditions. Because the model has to calculate overland flow and infiltration rates at the land surface, it was necessary to apply very short timesteps. The timestep size was adjusted automatically, and varied from a high value of 12 hours to as low as sixty seconds. To give an indication of computation time requirements, it took one hour of CPU time to simulate a six-year period on a Pentium 100 MHz computer.
5. Model calibration The model was calibrated over the period 1975 to 1995. This period includes both dry and wet years. No major flood impacts occurred, thereby providing a suitable period for the model calibration. Historical measurements of groundwater levels in more than 1200 bores screened in the sandy layers of the Shepparton Formation were used. Some data on streamflows and soil moisture was available for sufficient calibration of the channel and soil moisture components of the model. By including these processes in the model, we were able to find out if they play a major role in the water movement, and also to establish which data should be collected. The main emphasis of the calibration process has been to provide good temporal and spatial agreements between simulated and observed groundwater levels throughout the irrigation district. The calibration performance can be measured by the ability of the model to predict historical regional groundwater flow patterns and historic trends in groundwater levels over the calibration period. The results from the calibration runs are described in the following Sections 5.1–5.3.
Fig. 4. Comparison of observed and simulated groundwater levels at selected bores. The bore numbers are as they are stored in the database. They are located in Fig. 5 at coordinates (x,y,layer) as used in the model.
bore 287 over the period of 1975 to 1980 do not match the observed values. This can be attributed to flooding which, due to lack of information, was not taken into account in this model.
5.1. Graphs of groundwater levels vs time
5.2. Contour plots at selected times
The observed readings at different bores are compared with the ones computed by the model. Throughout the calibration process, a representative sample of 21 bores was selected for which the water levels were compared. Fig. 4 shows the time series of six such bores. The model results are shown as solid lines while the water level observations are shown as dots at the time the reading was taken. Fig. 4 shows that there is a general drop of the groundwater elevation which can be attributed to the groundwater pumping from the 48 pumps which have been operating since 1981. This observation is more obvious for the bores 230 and 389, which are situated within the area affected by the pumping scheme. From Fig. 4, it can be seen that the computed heads for
Fig. 5 presents a contour plot for May 1993. In Fig. 5 (a) the entire Wakool model is shown, while in Fig. 5 (b) a sub region of the model (Tullakool irrigation area) is presented in more detail. In Fig. 5, the evaporation basins are presented as well as the drainage channels. In order to evaluate the performance of the calibration of the model, plots similar to the one shown in Fig. 5 were plotted at approximately 5-year intervals. Since errors are accumulating over the calibration period, results at or near the end of the calibration period are expected to have maximum errors. For this reason, results at July 1994 (which is close to the end of the period of calibration) where selected to be presented in this paper. A calibration with an overall average error (difference
C. Demetriou, J.F. Punthakey / Environmental Modelling & Software 14 (1999) 129–140
135
Fig. 5. Comparison of observed (circles) and simulated groundwater heads at July 1994. (a) Model area, (b) Subsection of the model area showing in detail the evaporation basins (light blue), drainage channels (dark blue lines), extraction pumps (light blue dots). Colours show the groundwater levels as obtained from the model. The small circles indicated the observation bores, which were monitored at that particular time. A good calibration is achieved when there is a colour match between the colour in the circles and the colour of the surrounding area.
136
C. Demetriou, J.F. Punthakey / Environmental Modelling & Software 14 (1999) 129–140
between the computed results and observed values at any particular time at the observed points) of less than 1.0 m was achieved which can be considered as good taking into consideration the grid size of the model. Inspection of maps with simulated and observed water levels revealed that the model is able to achieve good agreement with the observed groundwater levels both spatially and temporally. It should be acknowledged that the developed model is based on various assumptions: boundary conditions; parameters for flow in the saturated and unsaturated zones; pumping stresses in the deeper aquifers; and others. The results for different scenarios produced by the model depend on these assumptions. Therefore, when more detailed information becomes available, the model should be recalibrated and the scenarios re-evaluated. 5.3. Percentage of the area affected by shallow watertable Another indicator, which was used in the model calibration, is the percentage of the area with watertable lying within 2.5 m from the ground surface. A comparison of the observed and computed percentage of the area affected by shallow watertable was carried out at July 1987 and 1995. The year 1987 is the year when the Wakool Sub-Surface Drainage (WSSD) scheme was fully operational and 1995 is the year when the Wakool LWMP commenced. Table 2 shows the computed results for these two landmark years as well as the corresponding observed values. As shown in Table 2, a high degree of agreement between the computed and the observed values has been achieved. 5.4. Sensitivity analysis of the Wakool model Sensitivity analysis was carried out to establish which parameters were most affecting the model. Model runs were carried out to investigate the model sensitivity on hydraulic conductivity, storage coefficient, pumping, leakage, and evapotranspiration.
Table 2 Comparison of the computed and observed percentage of the area affected by shallow watertable for July 1987 and 1995 % of area* in July 1987 Derivation of with watertable depths values (m BGL**)
Computed Observed *
% of area* in July 1995 with watertable depths (m BGL**)
⬍ 1.5
⬍ 2.0
⬍ 2.5
⬍ 1.5
⬍ 2.0
⬍ 2.5
9 9
13 12
18 16
5 6
7 7
12 13
The total model area is 810 (2 ⫻ 2 km) cells ⫽ 3240 km2; BGL: Below ground level.
**
5.4.1. Hydraulic conductivity The hydraulic conductivity was increased (from the spatially distributed calibrated values) by two orders of magnitude and decreased by two orders of magnitude but the change was in the order of a few centimeters up or down respectively which indicated that the effect was very small. This indicates that the process, which has the most pronounced effect on the watertable, is the recharge through the unsaturated zone. If the recharge through the unsaturated zone remains unchanged, the water level will not rise even with significant decrease in hydraulic conductivity. 5.4.2. Storage coefficient The storage coefficient for the confined layers was changed to values ranging from 0.5 to 1.25 times the calibrated spatially distributed values. The effect on the water level was insignificant (a few centimeters up or down) which indicates that the effect is minor. This again can be attributed to the fact that the recharge through the unsaturated zone is the process most affecting the watertable. 5.4.3. Pumping and channel leakage The pumping and channel leakage, were changed from 80% to 120% of their spatially distributed calibrated values. The effect on the water level was significant (up to 0.5 m change in general, and in some cases up to 1.0 m change) which indicates that the effect is locally significant. 5.4.4. Evapotranspiration The variation in evapotranspiration between 85% to 100% of the time-varying pan value had a significant impact on water levels (more than 1.0 m change). This indicates that this was a significant parameter. Field measurements of this parameter are needed in order to reduce the uncertainty in the calibrated model. Moreover, on-farm management options, which impact on evapotranspiration, may result in a significant change in the water levels. Based on the above descriptions, it can be said that the model is very sensitive to pumping, leakage and evapotranspiration. Pumping and leakage have a local effect while evapotranspiration has a regional impact.
6. Predictive simulations To determine the impact of changing management options over time, the timeframe for the model was extended by 25 years from 1995 to 2020. Before the commencement of the simulations, records (time series) on rainfall/irrigation, evapotranspiration, crop distribution, leaf area index, root depth and pumping/leakage need to be provided to the model up to the year 2020.
C. Demetriou, J.F. Punthakey / Environmental Modelling & Software 14 (1999) 129–140
These records are available up to 1995. Thus, artificial records were constructed to the year 2020. (a) Rainfall. The record was developed utilising the historical daily rainfall information at the Deniliquin meteorological station (70 km away from Wakool) extending from 1878 to 1996. A continuous 25-year portion of the record, which represented the average rainfall of the whole record, was established by utilising the moving average method with a moving average of 25 years. For an average rainfall over a 25-year period, the position of the average value in the moving average series indicates the relative position of the middle of the 25-year record in the actual rainfall record. This record was added on top of 1995 to extend the record to year 2020. A similar approach can be adopted for selecting a set of 25 “dry” years (by selecting the minimum value of the moving-average, or a set of 25 “wet” years (by selecting the maximum value of the moving average). (b) Irrigation. Based on the rainfall record and crop water demands, irrigation time series for each crop were developed. The crops used in the model are: no crop; rice; perennial (summer); and autumn crops. (c) Evapotranspiration. The evaporation data was obtained from the Deniliquin meteorological station, which has records from 1977 onwards. Thus, the evapotranspiration series for the selected 25-year period was synthesised using statistical methods based on the records of rainfall and evapotranspiration from 1978 to 1995. Knowing the evapotranspiration record, an actual evapotranspiration series for each crop was generated using the crop coefficients. (d) Crop distribution. The crop distribution of 1993/94 was used from 1993 onwards. The crop distribution was obtained from satellite images and consisted of the crops mentioned above. (e) Leaf area index and root depth. Both records were extended by repeating the existing yearly variations to year 2020. (f) Pumping/leakage records. In order to satisfy all the requirements by the different scenarios examined, various records for pumping and leakage rates were constructed. 6.1. Simulation of various management scenarios Various combinations of the following management options were simulated with the model. 6.1.1. Recycling ponds in conjunction with land forming (laser levelling) On-farm recycling ponds are specially constructed ponds in every farm, which will collect any surface runoff from rainfall and excess irrigation. For our modelling
137
exercise, the effect of the recycling ponds is simulated as recharge to the watertable caused by leakage. Due to the big grid size of the model, the effect of the recycling ponds is implemented in the model as points, which leak into the aquifer. The recharge rate was estimated to be 46 mm/year (Demetriou et al., 1998) and the area affected is proportional to the area of the farm. Laser levelling will be carried out in conjunction with the recycling ponds option. The aim of laser levelling is to contour the farm in such a way so that seepage will be minimised. Land forming (laser levelling) was implemented in the model by increasing the Manning surface number (Ritzema, 1994) in the model (the soil surface becomes smoother and therefore will be less seepage). The only information available on the implementation, was that the percentage of laser levelled area within the model, will start in 1995 and reach 45% at year 2000 and 75% at year 2005. Based on this, the Manning surface number was increased from 1 to 20 for randomly selected model grid cells in these two model periods. In addition, the detention storage coefficient (Ritzema, 1994) was modified. For the calibration period 1975 to 1995, it was set to 1.0 m for the cells covering the evaporation basins (see Fig. 2), 0.2 m for the rice cells and 0.01 m for the rest of the model area. From 1995 onwards, the detention storage coefficient for model cells other than those representing rice and evaporation basins, was decreased to 0.0075 m. 6.1.2. Deep rooted perennials The objective of planting deep-rooted perennials is to extract water through the root zone and release it to the atmosphere by transpiration. The effectiveness of this scheme is limited by the root-depth and the salt tolerance of the plant. For this model scenario, the root depth of the summer (perennial) crops was increased from 0.5 m to 2.0 m. This option becomes active in the model from 1995 onwards. 6.1.3. Tree planting The objective of tree planting is to capture, through the root zone, any freshwater which might leak from irrigation channels. The effectiveness of this scheme is depended on the density and depth of the roots as well as the salt tolerance of the trees. Line extractions were placed along the irrigation channels (100 m on either side) within the Tullakool Irrigation area operating at an extraction rate of 1.0 mm/day (200 m3/day/km—strip). The width of the strip is 200 m, 100 m on either side from the irrigation channels. The water extraction was dropped to 35% during the winter months. The trees are modeled as pumps, and they extract water over the whole depth of the top layer (Upper Shepparton Formation). Tree planting becomes active from 1995 onwards and the trees are assumed fully mature from the time of planting.
138
C. Demetriou, J.F. Punthakey / Environmental Modelling & Software 14 (1999) 129–140
6.1.4. Shallow pumping The objective of shallow groundwater pumping is to pump the salty groundwater to the evaporation basins and thus slow or reverse the rise of the watertable. An extra, 48 shallow pumps were added to the model and were installed within the Tullakool irrigation area. Their locations were nominated by the office of DLWC in Deniliquin. The pumping rate for each pump is 250 ⫻ 103 m3/year starting from 1995. 6.1.5. Deep groundwater pumping options The objective of deep groundwater pumping is to reverse the upward movement of salty water from the lower to the upper aquifer. The effect of deep groundwater pumping is expected to be a long term one. Three deep pumps (extract water from the Calivil/Renmark Formations) were installed within the Tullakool irrigation area at the locations nominated by the office of DLWC Deniliquin office. The rate is 7.3 ⫻ 106 m3/year and starts from 1995 onwards. 6.2. Results of scenario modelling Table 3 presents a summary of the results from the different scenarios examined. The percentage of the area with watertable lying within 1.5 m, 2.0 m and 2.5 m from the land surface is presented in Table 3. Scenario 1 shows the results for the “basic” (as in 1995) scenario for the years 2005 and 2020. At year 2020, the area with watertable within 2 m reaches 29% of the total area. Scenarios 2 to 6 present the results obtained from recycling pond in conjunction with other options. Under these scenarios, the best performance (the least percentage of the area with the watertable lying within
2.0 m of the ground surface) was achieved by scenario 6 (48-pump option). Scenarios 7 and 8 cover the cases where all the options are incorporated at different implementation rates and with different detention storage. In scenario 7 the detention storage coefficient reduced to 50% of the ‘as in 1995’ case, while in scenario 8 it was reduced to 75%. As expected, slightly better performance was achieved by the smaller detention storage. Scenarios 9 and 10 investigate the effect the rainfall variation has on the results. The results are very sensitive to the average annual rainfall but more importantly to rainfall intensity and distribution with respect to time. Scenarios 11 to 14 investigate the effect of the different options as if they were implemented individually. The best performance of these was obtained from the 48-pumps option (scenario 13). 7. Model audit by external reviewers External consultants reviewed the developed model and a report with comments and recommendations from the reviewers as well as from the modelling team, was produced (Meynink et al., 1997). The following are some of the recommendations from the review process: 쐌 Field measurements need to be carried out in order to establish better estimates for the important parameters. These parameters are: evapotranspiration and leaf area index; surface roughness coefficient; detention storage coefficient. For deep-rooted perennials, information is needed on root development, root density and tolerance to groundwater salinity, and for tree planting, capacity of the plant to extract water, rates of growth of tree and salt tolerance.
Table 3 Percentage of the area affected by shallow watertable under various simulated management scenarios for March 2005 and 2020 Simulated scenario
% of area* in March 2005 with watertable depths (m BGL**) ⬍ 1.5
1. Basic (as in 1995) 2. Rec. Ponds ⫹ (land forming) LF1 3. Rec. Ponds ⫹ (land forming) LF2 4. Rec. Ponds ⫹ LF1 ⫹ deep roots 5. Rec. Ponds ⫹ LF1 ⫹ trees 6. Rec. Ponds ⫹ LF1 ⫹ 48 pumps 7. Rec. Ponds ⫹ LF1 ⫹ 48 pumps ⫹ 3 deep pumps ⫹ deep roots 8. Rec. Ponds ⫹ LF2 ⫹ 48 pumps ⫹ 3 deep pumps ⫹ deep roots 9. Rec. Ponds ⫹ LF2 ⫹ 48 pumps ⫹ 3 deep pumps ⫹ deep roots (dry) 10. Rec. Ponds ⫹ LF2 ⫹ 48 pumps ⫹ 3 deep pumps ⫹ deep roots (wet) 11. Basic (as in 1995) ⫹ deep roots 12. Basic (as in 1995) ⫹ trees 13. Basic (as in 1995) ⫹ 48 pumps 14. Basic (as in 1995) ⫹ 3 deep pumps
5 3 4 3 2 2 2 2 1 3 5 3 3 3
⬍ 2.0 9 6 7 6 5 3 2 2 2 7 7 5 4 5
⬍ 2.5 12 10 11 10 9 7 5 5 4 16 12 10 9 10
% of area* in March 2020 with watertable depths (m BGL**) ⬍ 1.5 20 7 10 7 4 2 2 2 1 2 13 9 1 10
⬍ 2.0 29 10 14 10 8 4 3 3 1 4 23 17 17 19
⬍ 2.5 40 16 22 14 12 8 3 8 2 10 34 23 27 28
The total model area is 810 (2 ⫻ 2 km) cells ⫽ 3240 km2; LF1: Land forming with detention storage coefficient reduced by 50% of the "as in 1995" case; LF2: Land forming with detention storage coefficient reduced by 75% of the "as in 1995" case; **BGL: Below ground level.
*
C. Demetriou, J.F. Punthakey / Environmental Modelling & Software 14 (1999) 129–140
쐌 Continue water level and water quality monitoring, preferably at quarterly intervals, or, at least take readings twice a year. 쐌 Expand the deep bore network and continue water level and water quality monitoring (minimum twice a year, preferably quarterly). 쐌 Field verifications of the leakage locations and rates need to be conducted. 쐌 In order to obtain an indication on the validity of model outputs, uncertainty analysis should be carried out. This requires a great amount of extra work because many different models need to be set up and calibrated.
139
From the results obtained, it can be concluded that the model is very sensitive to parameters such as evapotranspiration, rainfall pattern, and the parameters associated with land forming (Manning surface number and the detention storage). The results are strongly influenced by the above parameters and field measurements need to be carried out in order to establish better estimates for them. In addition, the rise and fall of the watertable is strongly influenced by the intensity, distribution and duration of the rainfall. Despite the data limitations identified in Section 7, the use of MIKE SHE to develop a model for the Wakool Irrigation District has provided an important tool for evaluating possible future strategies for watertable reduction in this area.
8. Conclusions and recommendations The MIKE SHE modelling package was chosen over other packages such as the widely used MODFLOW, because it is well suited to describe the dynamic interaction between the surface and subsurface water systems, which is particularly important in the high watertable areas of the Wakool Irrigation District. A range of groundwater management options proposed in the land and water management plan (LWMP) for the Wakool Irrigation District have been investigated and their effectiveness have been evaluated. These options range from revegetation-type alternatives such as planting trees, salt tolerant deep rooted perennials, to engineering alternatives such as pumping from different aquifers, laser levelling (land forming), and utilising recycling ponds. The model predicts that with the current irrigation practices, the area with a watertable less than 2 m deep will further expand, and by the year 2020, will cover about 29% of the Wakool Irrigation District, an increase of 22% compared to that at 1995. Of all the groundwater management options assessed in this study, on-farm recycling in conjunction with land levelling (scenario 3) and the 48-pump option (scenario 13) perform the best. Under these two scenarios, the shallow watertable, at year 2020, is expected to reach 14% and 17% of the area respectively (an increase of 7% and 10% compared to 1995). It should be pointed out that the confidence in the 48-pump option (scenario 13) is very high because the pumping rates are known. On the other hand, the results from the on-farm recycling option (scenarios 2 and 3) carry a high degree of uncertainty since the coefficient for surface roughness and detention storage are not well known. Also, high uncertainty is associated with the deep-rooted perennials option (scenarios 4 and 12). This is due to many factors associated with the root development such as root density, root depth and plant tolerance to groundwater salinity.
Acknowledgements Mr. Borge Storm, consultant, Lawson and Treloar Pty Ltd, is acknowledged for his input in setting up the model and providing ideas on implementing and analysing the scenario simulations. Mr. K.R. Bogoda and Mr. N. Kulatunga, DLWC, Deniliquin are thanked for their advice and provision of field data. Mr. H. Schroo, DLWC, Parramatta is acknowledged for his constructive criticism and input during the model development and analysis of the scenario simulations.
References Bogoda, K.R., Kulatunga, N., Hehir, K., 1994. Overview of hydrogeology and assessment of subsurface drainage option for watertable control in Wakool. Internal Working Document, DLWC, pp. 123. Carr, R.S., Punthakey, J.F., Cooke, R., Storm, B., 1993. Large-scale catchment simulation using the MIKE-SHE model: 1. Process simulation of an irrigation district. International Conference on Environmental Management. Wollongong, Australia, Balkema, Rotterdam, pp. 467–472. Danish Hydraulic Institute, 1993. MIKE SHE water movement-user’s guide and technical manual, ed. 1.0. DHI, Denmark, pp. 81. Demetriou, C., Storm, B., Punthakey, J.F., 1998. Development of an integrated model to simulate on-farm water recycling in the Wakool Irrigation District. Second International Conference on Environmental Management (ICEM2), vol. 1. Wollongong, Australia, pp. 633–640. Kristensen, K.J., Jensen, S.E., 1975. A model for estimating actual evapotranspiration from potential evapotranspiration. Nordic Hydrology 6, 70–88. McDonald, M.G., Harbaugh, A.W., 1988. A modular three-dimensional finite-difference ground-water flow model. U.S. Geological Survey techniques of water resources investigations, Book 6, pp. 586. Meynink, W.J.C., Doherty, J., Demetriou, C., Punthakey, J.F., 1997. NSW DLWC Wakool SHE model review. Report # psm— wakool—051.1.doc, Pells Sullivan Meynink (Brisbane). Ritzema, H.P., 1994. Drainage principles and applications. International Institute for Land Reclamation and Improvement (ILRI) publication 16, 2nd ed. The Netherlands, pp. 1125.
140
C. Demetriou, J.F. Punthakey / Environmental Modelling & Software 14 (1999) 129–140
Storm, B., Punthakey, J.F., 1995. Modelling Environmental changes in the Wakool Irrigation District. MODSIM 95, Int. Congress on
Modelling and Simulation Proceedings, vol. 3. Newcastle, Australia, pp. 47–52.