Received 14 January 2016; accepted 14 February 2016; published 17 February 2016

1. Introduction

Management of regional groundwater systems requires an understanding of the balance between recharge, inter- aquifer leakage and lateral movement of groundwater. There are currently few reliable hydraulic or tracer based methods for estimating rates of diffuse inter-aquifer leakage through an aquitard. The tracer-based approaches involve one or more of the following: environmental or injected tracer techniques [1] ; advection-diffusion modelling of pore water chloride profiles [2] -[4] , adoption of the compartmental mixing model approach [5] . Hydraulic based approach to the assessment of site specific vertical leakage from Darcy’s equation involve: vertically averaged conductance [6] ; based on vertical hydraulic gradient and Darcy flow equation [7] ; equivalent hydraulic conductivity approach [8] ; adaptation of radial flow theory from pumping test data [9] ; analytical solution to pressure perturbation and fluid leakage through the aquitard [10] . Hydraulic-based quantification of vertical leakage can be difficult since the key parameter, vertical hydraulic conductivity can vary considerably. Similarly, tracer or hydraulic based site specific recharge or vertical leakage estimation becomes problematic unless duality of flow regimes (diffuse and preferential) is considered [11] [12] .

In the Glencoe-Nangwarry-Nagwarry (GNN) recharge zone (Figure 1) of the Gambier Basin, South Australia, Brown et al. [13] used environmental isotope techniques to study the dynamic behavior of the groundwater systems and concluded that vertical leakage to the confined aquifer is attributed to preferential flows (fractures, faults or sinkholes) rather than via porous media flow through the aquitard. Recently, Somaratne et al. [14] extended the study by Brown et al. [13] to an additional 6 sites and have shown that vertical leakage to the confined aquifer in the GNN recharge zone is characterized by both diffuse and preferential flow processes.

Transient numerical groundwater models are widely used to provide a regional or sub-regional estimate of vertical leakage through aquitards [15] by inverse calibration to heads and groundwater flow rates [16] as models are independent of the mode of recharge and leakage processes. Several groundwater models have been developed for the entire Gambier Basin [17] [18] or outside the GNN recharge zone [19] . These models are at their conceptual stages and may require further work to be used as tools for regional or local water balance calculations.

Groundwater models are valuable tools for analyzing groundwater systems. Nevertheless, applying groundwater models to field problems can result in errors from conceptual deficiencies, numerical errors, and inadequate parameter estimation. As such, models cannot be validated [20] and viewed as the application of distributed hydrological models is more an exercise in prophecy than prediction [21] . Thus, the primary value of models is heuristic [22] . In general, manual procedures for groundwater model calibration can be extremely time-consuming and frustrating [23] , and inverse modelling is a key step in groundwater related studies [24] . Boyle et al. [25] however, suggests using the combined strengths of the manual calibration approach, which has been developed and refined over the years; and automatic calibration. The study of Saiers and Bolster [26] shows that even relatively simple models can predict transient groundwater flow in a complex aquifer with reasonable accuracy. Model calibration allows reduced parameter uncertainty and, therefore uncertainty in simulation results [27] . According to Moore and Doherty [28] , no matter which regularization methodology is employed for model calibration, the inevitable consequence of its use is a loss of detail in the calibrated field leading

to erroneous predictions made by a model that is ostensibly “well calibrated”. Yeh et al. [29] recently addressed the uniqueness, scale, and resolution issues in groundwater flow models. They emphasized the need for full specifications of; flux boundaries, sources and sinks, and heads everywhere in the model domain, to obtain a unique solution to transient flow. The impracticality of such a requirement confirms the heuristic suggestion of Oreskes et al. [22] .

Most recharge estimation efforts within the Gambier Basin [30] -[34] (Figure 1) have involved the shallow unconfined Tertiary limestone aquifer (TLA) as the TLA is the main source of irrigation, industrial, livestock and town water supplies. The Tertiary confined sand aquifer (TCSA) in the Gambier Basin is also an important groundwater resource that supplies good quality water to both urban and rural users, including water supplies to seven townships. Licensed groundwater extraction for irrigation in 2013-14 from the TCSA was 16 × 10⁶ m³ [35] . Since the confined aquifers are sensitive to changes in flow rates that may result in excessive drawdowns, extraction practices or reduction in vertical leakage in recharge zones can have a detrimental effect on groundwater levels and water quality. To achieve sustainable groundwater management, it is necessary to understand how, and by how much, TCSA receives vertical leakage via diffuse and preferential flows from the overlying TLA through the aquitard.

Recharge to the TLA and vertical leakage to the TCSA has spatial and temporal variations. Relatively little work has been done on vertical leakage to the TCSA and no quantification of recharge in the GNN area has been completed. The present study adopts a heuristic approach for sub-regional estimate of vertical leakage through the aquitard. A groundwater flow model spanning the Glencoe and Nangwarry areas of the GNN recharge zone was developed to improve understanding of recharge to the TLA and vertical leakage to the TCSA aquifers and to assess the impact on the regional flow system.

2. Hydrogeological Setting of the Model Domain

Since a detailed description of the Gambier Basin and characteristics of the GNN recharge zone is given in Somaratne et al. [14] , only a summary is provided here, focusing on features relevant to vertical leakage processes. The model domain spans parts of the Glencoe, Nangwarry and Penola areas (Figure 1). Numerous structural features occur within the model domain in terms of orientation and distribution of faults which are important elements for vertical leakage. Two major aquifer systems considered within the model are the TLA and the TCSA. Karstic features are common within the Gambier limestone, and it is a dual porosity system. Lithologies of the upper Tertiary aquitard overlying the TCSA are glauconitic and fossiliferous marls and clays of the Narawarturk Marl and clays of the Mepunga Formation [36] . The TCSA is a multi-aquifer system but due to the paucity of data, is treated as one aquifer unit. Since most of the wells drilled penetrate only the upper sand unit, hydraulic interconnection between the sub-aquifers is unknown [19] .

The area is highly faulted. Groundwater bulging occurs in the TCSA within the GNN recharge zone, with dominant mounds at Glencore, SA 4 and VIC 3 sites (Figure 1), indicating that TCSA receives high vertical leakage. In this area, Love [36] , Brown et al. [13] and Somaratne et al. [14] reported the presence of high percentages of radiocarbon activity indicating the TCSA receiving modern recharge. About 83% of the Hundred of Nangwarry area is covered with plantation forest, mainly Pinusradiata. In general, groundwater movement from the GNN recharge zone spreads outwardly in a northwesterly to southerly direction. The topography of the area is generally flat at about 70 m above Australian Height Datum (m AHD). The main geological units in a downward order consist of the unconfined TLA, Narawarturk Marl aquitard, Mepunga Formation, Dilwyn Clay aquitard and the TCSA (Figure 2).

Somaratne et al. [14] described eight investigation sites in this area; four (SA 1, SA 2, SA 3, SA 4) are located within the South Australian side of the GNN recharge zone and the rest (VIC 1, VIC 2, VIC 3, VIC 4) are located across the State boundary on the Victorian side (Parish of Nagwarry) of the recharge zone (Figure 1). The main features of the study sites include: the presence of Narawarturk Marl; 9 m thick at SA 3, 1 m thick at VIC 1, and at VIC 4 about 12 m thick; varying thickness of Mepunga Formation (4 - 12 m at SA sites and 1 - 12 m at VIC sites) and Dilwyn Clay aquitard (2 - 7 m at SA sites and 3 - 20 m at VIC sites). Overall, no aquitard occurs at VIC 3 (only 3 m Mepunga Formation separates TLA and TCSA) and while it is thickest at the VIC 4 site (26 m: 9 m Narawarturk Marl, 12 m Mepunga Formation and 5 m Dilwyn Clay). Apart from the SA 1, VIC 3 and VIC 4 sites, all sites are located within plantation forest.

The presence or absence of a competent aquitard and the head difference between the TLA and the TCSA determines the vertical leakage rate and is an important consideration in model conceptualization. In the south- west part of the GNN recharge zone, the Narawarturk Marl and Mepunga Formation are thin. For example at HIN079 both the Narawarturk Marl and Mepunga Formation are about 1 m thick, and the 12 m thick Dilwyn Clay aquitard is present. About 20 m uplift of all Formations is evident from drillhole 7022 - 8033 through to

drillhole 7022 - 6448. A major displacement occurs between two fault lines east and west of TWS 3, the water supply well to the township of Kalangadoo. The thickest TLA is found at this location with thin Narawarturk Marl and Mepunga Formation. The area is identified as the Kalangadoo-Tarpeena fault zone. In this area the Narawarturk Marl and Mepunga Formations are about 1 m thick, while the Dilwyn Clay aquitard is about 20 m thick. Vertical displacement of the Formations combined with the thin aquitard favours relatively high direct recharge to the TCSA, as is evident in the presence of the groundwater mound in the Glencoe area (Figure 1). Towards the east of the Kalangadoo-Tarpeena fault zone, from NAN042 to NAN084 wells, a moderately thick Narawarturk Marl (3 - 7 m) and Mepunga Formation (3 - 12 m) are found. Monitoring well NAN042 is the TCSA well at study site SA 4 and is located on a fault line. Monitoring wells TLY011 and TLY006 are the TCSA wells at study sites VIC 2 and VIC 1. In this area once again, the thickness of the Narawarturk Marl (1 - 2 m) and Mepunga Formation (3 - 5 m) is reduced, but thick Dilwyn Clay (13 - 20 m) is present. Investigation sites VIC 1 - 4 are outside the model domain.

3. Methods

The results from Somaratne et al. [14] show that TCSA receives both diffuse and preferential vertical leakage from the GNN recharge zone. A transient groundwater model previously developed in 2007 for the Nangwarry and Penola wellfield protection zone delineation is used to estimate recharge to the TLA and vertical leakage to the TCSA in the South Australian portion of the recharge zone. The model setup comprised two layers corresponding to TLA and TCSA with 100 m × 100 m grid spacing, with 300 rows and 400 columns. The model is based on MODFLOW-2000 [37] , and has a simulation period of 20 years from October 1985 to September 2005. Spatially interpolated cell-by-cell, initial head distributions were obtained by kriging of observed head values of monitoring wells. The model uses monthly stress periods to simulate temporal behaviour of well hydrographs. A quasi-three-dimensional approach [38] is adopted, where resistance to flow in low conductivity units (Narawarturk Marl, Mepunga Formation and Dilwyn Aquitard) is lumped together into a vertical leakance term between the adjacent layers and is represented by vertical flow only with no storage effects. Here, Mepunga Formation is a low yielding minor aquifer which has not been developed within the Gambier Basin, except at Kingston (Figure 1) town water supply. Horizontal flow may exist in the Mepunga Formation within the model domain, but is considered negligible compared to flows within the TLA and TCSA.

General head boundary (GHB) conditions of varying conductance were used around the model domain to represent the external hydrological conditions. External heads forcing flux into or out of the model domain at GHB boundaries were set according to time averaged groundwater levels outside the model domain. Conductance parameters along the GHB boundaries varied according to the updated hydraulic conductivity along boundary cells as the calibration progressed. Within the model area, hydrological stresses such as recharge, evapo transpiration and groundwater extraction were used. Layer elevations were defined using lithological information from drillholes within and surrounding area of the model domain; 120 drillholes were used for unconfined aquifer elevations, 108 for aquitard and 70 for TCSA. Once the layer elevations were identified, the next step in model setup was to define aquifer parameters; horizontal hydraulic conductivity (K_h) for the TLA was obtained from three-pumping tests and 11 estimates based on specific capacity-transmissivity relation established by Mustafa and Lawson [39] . For TSCA, K_h was available for nine locations from 72 hour pumping tests, and specific storage (S_s) for two locations. The anisotropy for K_h to K_v (vertical hydraulic conductivity) of aquifers were taken as 10:1, for TLA and TCSA. For K_v of the Dilwyn Clay aquitard hydraulic conductivity was estimated using particle size distribution data using the quasi-numerical model, ROSETTAV1 [40] . The core samples for analyses were taken from Nangwarry area, and show relatively high percentage of silt and sand (% Clay = 32 - 42, % of Silt = 20 - 53, % of Sand = 15 - 48). For comparison, initial estimates of parameters, calibrated parameters and parameters derived from post-modelling pumping test conducted in 2011 within the recharge zone [41] are given in Table 1.

Apart from the town water supply wells, actual groundwater extractions volumes for irrigation were not available and therefore, licensed allocations of the irrigated areas were used. A total of 28 × 10⁶ m³∙year⁻¹ crop area-water requirement (based annual allocations) were used as the groundwater extraction from the TLA aquifer at licensed property locations. About 26% of the model domain is covered with plantation forests; pine (Pinusradiata) and Tasmanian blue-gum (Eucalyptus globulus). Evapotranspiration (ET) is an important sink for shallow groundwater systems, and simulates the effects of plant transpiration, direct evaporation and removing water from the saturated groundwater regime. For shallow groundwater, evapotranspiration (ET) was repre- sented using the EVT package of the MODFLOW. Tree and soil-based ET measurements have been undertaken by the former CSIRO Division of Forestry (Personnel communication with Richard Benyon, CSIRO) and ET data for the study was obtained for the Nangwarry plantation forests (Figure 3). For that part of the plantation forest burnt during the 1983 bush-fire, the year of canopy closure was taken as 1991. Effective rooting depth for plantation forest was set as 6 m [42] .

For effective model calibration, it is essential to capture the spatial and temporal variation of recharge in the study area. An initial map of long-term recharge distribution was prepared using recharge values obtained by the water table fluctuation (WTF) method for 42 locations within and outside the study area. Recharge values ranged from 4 - 216 mm∙year⁻¹. The model area was divided into 11 recharge zones. Following Hanson et al. [43] , a rainfall-recharge relationship was established with recharge as a percentage of infiltration of excess rainfall to the groundwater system. In the area, daily maximum ET for winter months is about 2 mm and for summer is about 5 mm. In the study area, approximately 25 mm (6%) of monthly rainfall comprised of daily rainfall < 2 mm during winter months and that for summer months, about 41 mm (23%) of monthly rainfall comprised of

^*From pumping tests outside the model domain; Post modelling data in 2011 for TCSA recharge zone at SA 1 - 4 sites.

daily rainfall < 5 mm. It is considered that daily rainfall less than daily ET would not contribute to soil infiltration. These rainfall threshold limits were deducted from respective winter and summer monthly rainfall to obtain monthly excess rainfall, and a fraction of that would potentially contribute to recharge.

Therefore, a monthly recharge estimate is obtained from:

(1)

where RCH is the monthly recharge, R_Month is the monthly rainfall, R_Threshold is the threshold rainfall and, x is the fraction (%) that would recharge. The percentage (x) is determined through trial-and-error model calibration for different recharge zones. The calibrated x values range from 5% for Penola forest to 85% at Glencoe groundwater mound area, maintaining the recharge range obtained using the WTF method.

One of the important issues considered in model calibration is to limit the number of parameters to be estimated without losing too much of the spatial variability of parameters derived from pumping tests and recharge estimates from the WTF method. The exception to this is groundwater mounds locations in the TCSA and locations of major faults. For parameters and recharge, the well-established zonation approach was adopted [44] [45] . The model was calibrated to represent observed areal head distributions and the well hydrograph response due to natural and seasonal variation, ET and pumping stresses at specific sites. The calibration procedure involved modifying specific yield, specific storage, vertical leakance and horizontal hydraulic conductivity to improve regional hydraulic gradient and potentiometric contour patterns.

4. Results and Discussion

4.1. Transient Model Calibration

The calibrated groundwater potentiometric heads for the end stress period were compared with the measured heads and are shown in Figure 4(a) and Figure 4(b).

Model performance can be examined more rigorously by analysing the hydrographic responses at specific sites. At this stage aerially calibrated hydraulic conductivity and specific yields/storage were used as the initial parameter values. Local changes were made to parameters to reflect the hydraulic gradients and observed bore hydrographs. Following this, once again, parameter values were re-calibrated to match potentiometric contours. This iterative process continued until satisfactory areal and bore hydrograph calibration was achieved. Calibrated bore hydrographs within plantation forests and irrigation areas are depicted in Figure 5 and Figure 6.

We calibrated all long-term data available from monitoring wells in the model domain. Twelve monitoring wells were used within the plantation forest for TLA calibration. Apart from the PEN002 and PEN011 wells, which are influenced by groundwater extraction from adjacent irrigation areas, the observed seasonal range of water level fluctuations within the plantation forest were minor, and long-term observed and calculated water level decline agree reasonably well. In contrast, observed groundwater level fluctuations vary substantially in the irrigation area (Figure 6) as some of the pumping wells are being used as head observation wells (Figure 7), thus reducing the credibility of these observations.

When low credibility observations are present, difficulties arise as to overall model calibration. One problem is that time series do not display distinct seasonal and annual variations. This is further exacerbated due to water level monitoring undertaken twice per year, in summer and winter, yielding water level data with extreme values. Therefore, we attempted to calibrate the “average long-term trend” of the water table fluctuations of observed time series.

Paired TLA-TCSA wells, located close to each other, were calibrated to obtain reliable vertical leakage to the TCSA (Figure 8). Wells PEN008-PEN025 are located in the Penola forest; NAN055-NAN043 and NAN046- NAN042 are located in the Nangwarry forest. MON18-MON019 (Monbulla area), YOU029-YOU034 (Young area) and MIN016-MIN027 (Mingbool area) were available and used for irrigation areas. Apart from the over prediction at NAN042, head differences were reasonably well simulated through the modelling process.

Despite the groundwater mound in the Glencoe area spanning from Glencoe to Kalangadoo in the GNN recharge zone (Figure 1), no paired wells are currently available in the entire western region. Two TCSA wells in the Kalangadoo-Tarpeena fault zone were calibrated (Figure 9), but once again head-observations from these wells are of low credibility as GRY021 is equipped with a submersible pump and used for irrigation and livestock watering. The well GRY016 is currently backfilled and replaced, but was the water supply well for the Kalangadoo township during the study period. This Glencoe to Kalangadoo portion of the model heavily depended on areal calibration of head contours in the TCSA.

Quantitative calibration performance measures can be related to calculation of potentiometric head residuals and associated statistics. Figure 10 shows the scatter plot of calculated versus observed groundwater heads for TLA and TCSA. The total number of observations for the 20-year simulation is 2375.

Small positive skewness (0.2 m) close to zero is considered a reasonable calibration with a root mean squared error (RMSE) of 0.83 m. The normalized RMSE (expressed in %) is considered a more representative measure of fit since it accounts for the scale of the potential range of data values [46] , which is 2.2% in our case. In general, a normalized RMSE value less than 5% is considered a good fit between calculated and observed values. The correlation coefficient of the observed heads to calculated heads in this calibration is 0.995, indicating that the two range of data (calculated and observed) move together.

Sensitivity analyses were carried out to explore and quantify the impact of possible errors in input data on model performance indices. This also can be used to provide a first order uncertainty analysis [47] . Sensitivity analyses for the key parameters, K_h, K_v, S_y, S_s, and recharge were undertaken. Each set of values was altered separately by increasing and decreasing the calibrated parameter value using multiplication factors between 0.1 and 10. The parameter multiplier is spatially and temporally constant. The results of sensitivity analyses show hydraulic conductivity and recharge to TLA are most sensitive (Figure 11). The key hydraulic parameters determined during the post-modelling period in 2011 are in good agreement with the range of values obtained from model calibration (Table 1). The exception to this finding is the increased K_v of the aquitard at locations of faults and groundwater mounds in the Nangwarry and Glencoe areas.

4.2. Water Balance

Figure 12 shows the transient water balance for the TLA aquifer. Recharge to the aquifer varies significantly, 0 - 2677 m³∙day⁻¹, with an average 455 m³∙day⁻¹. The highest recharge occurred in response to June 2003 monthly rainfall of 184.6 mm, about 30% of the average annual rainfall in the area. Since recharge is a function of monthly rainfall this variation is expected. In response to the rise and fall of water levels, TLA storage changes with a net average depletion of 47 m³∙day⁻¹. Other flow components of the water balance include: leakage to TCSA (297 ± 11 m³∙day⁻¹), well extraction (79 m³∙day⁻¹), ET (36 m³∙day⁻¹) and net boundary flow out from the domain at 79 m³∙day⁻¹, all of which remained relatively stable over the study period.

Recharge to TLA and vertical leakage to TCSA was assessed on the basis of the water balance of management areas of the water allocation plan (WAP) for comparison with WAP recharge (Figure 13). Recharge for the current WAP was assessed using the WTF method [48] except for the Hundred of Young, which is based on previous studies [48] using environmental chloride and tritium. For simplicity, a specific yield value of 0.1 was taken for recharge estimation for the WAP.

WAP adopted recharge for Penola, Monbulla and Young are higher than the recharge from the calibrated model. For the Hundred of Grey, a single monitoring well was used to calculate a recharge at 150 mm∙year⁻¹ for the WAP [48] , while the model calculated recharge from areal calibration of 219 mm∙year⁻¹. Model calculated recharge (150 mm∙year⁻¹) for the irrigation areas of Young and Mingbool are equal since a single water balance zone was used in the model.

For the non-forested area of the Hundred of Penola, the model calculated an annual recharge of 87 mm∙year⁻¹, while for the Penola forest recharge was calculated at 19 mm∙year⁻¹ for TLA. When compared to non-forested areas, a recharge reduction of 78% is observed for the Penola forest. Despite the fact that 83% of the Hundred of Nangwarry is covered with plantation forest, a comparatively high recharge (80 mm∙year⁻¹) to TLA is calculated by the model. This represents 44% less recharge than adjacent non-forested areas within the Hundred of Nangwarry. Both Penola and Nangwarry areas are faulted (Figure 1).

The hydrogeological difference between Penola and the Nangwarry forest is the presence of preferential flow paths via a generally thin aquitard in Nangwarry. In contrast, Penola and Monbulla have about a 50 m thick aquitard layer (combined Narawarturk Marl, Mepunga Formation and Dilwyn Clay) separating the TCSA from the TLA. This indicates that the influence of structural faults is effective only in areas where the aquitard is thin or absent [14] . As such, the model calculates 8 mm∙year⁻¹ vertical leakage to TCSA in the Penola forested area, compared to 84.5 mm∙year⁻¹ vertical leakage in the Nangwarry forested area. One point to note is that the model results show leakage to TCSA as higher (84.5 mm∙year⁻¹) than the recharge to TLA (80 mm∙year⁻¹) in the Hundred of Nangwarry. This, results in a continuous decline of water levels in TLA. Combined with a groundwater loss to ET (13 mm∙year⁻¹) from areas where depth to water < 6 m, this would result in the TLA beginning to dry up. This phenomenon appears to be happening at investigation site SA 2 [14] .

Our results compare favourably with research indicating that groundwater recharge under Pinusradiata plantations is either zero [30] [31] [49] or substantially lower than under grassland [32] . Direct measurements of water use by plantations confirm recharge is usually zero after canopy closure [42] [50] [51] . However, the above cited studies either, do not consider preferential recharge to groundwater or used measurements made outside the GNN recharge zone. Compared to the average recharge from agricultural lands, Benyon and Doody [42] estimate that the reduction of recharge under Pinusradiata is about 85% and under Eucalyptus globlus about 78%. The model calculated a recharge reduction of 78% in Penola forest, where the recharge process is predominantly diffuse.

As noted by Bevan and Young [52] , forecasting future management conditions is heavily reliant on several assumptions, particularly those related to rainfall dependent recharge, groundwater use and other land use practices. Therefore, “what-if” scenarios or projections were considered to further assess the impact of land use on groundwater recharge. A simulation of replacing plantation forest with dryland pasture in the Hundred of Nang- warry area for the calibration period resulted in recharge to TLA of approximately 135 mm∙year⁻¹ and leakage to TCSA, of 98 mm∙year⁻¹.

With respect to replicating head contours and observed time series, the performance of our model may be considered appropriate for assessing vertical leakage to TCSA and the impact of the plantation forest on groundwater recharge. One of the major difficulties in model construction was to assign reliable K_v distributions, since the only data sets available at the time, were particle size distribution based estimated values (Table 1) and aquitard hydraulic conductivity (range 10⁻³ - 10⁻⁷ m∙day⁻¹) obtained from triaxial tests [36] . In addition, the ranges of transmissivity for TLA (200 - 10,000 m²∙day⁻¹) and TCSA (200 - 1600 m²∙day⁻¹) in the Gambier Basin [36] were used as a guide. Despite the calibration approach chosen, vertical leakage is truly vertical only if; the head in the layer supplying the leakage is constant, and if permeability contrasts in the aquitard and confined aquifer, and storage in the aquitard layer are negligible [53] . These assumptions may not fully hold in the GNN recharge zone. Harrington and Cook [54] highlight the needs for applying correction for measured confined aquifer pressure levels for the mechanical loading and unloading effect of the watertable levels in overlying unconfined aquifer. The study of Harrington and Cook [54] shows that the maximum possible loading/unloading impact on a confined aquifer can only be a fraction of the observed watertable fluctuation of the unconfined aquifer, equivalent to the specific yield. In this study, TCSA head corrections were not made and hence calculated pressure levels may include smaller over and under estimation than actually observed.

The use of groundwater allocation as actual extraction is another source of error. During the post modelling period, metered annual extraction data were available for the model area. Compared to the 2013-2014 irrigation season, actual groundwater extraction is 69% of the area-based annual allocation. This means an excess of 8.68 × 10⁶ m³ of water was extracted in the TLA water balance, which results in higher calculated recharge. Distribution of estimated excess volumes across the total irrigation area, 970 km², results in approximately 9 mm over estimation of recharge to TLA.

5. Conclusions

A groundwater flow model spanning the Glencoe and Nangwarry areas of the GNN recharge zone was developed to improve understanding of vertical recharge to the TLA and TCSA aquifers and to assess the impact on the regional flow system. The model was developed using MODFLOW-2000 with a calibration period from October-1985 to September-2005. Despite many studies considering the importance of vertical leakage to the confined aquifer, and the impact of plantation forestry on groundwater recharge, quantification of vertical leakage to the confined aquifer was lacking.

The model calculated recharge to the confined aquifer in diffuse recharge zone under plantation forest (Penola forest) of about 8 mm∙year⁻¹, a 78% reduction of recharge compared to adjacent non-forest areas. Where structural uplift and geological faults are present as in Hundred of Nangwarry, vertical leakage is about 84.5 mm∙year⁻¹. The higher recharge is attributed to the thin or absent aquitard and preferential flow through faults. The results gained in this study show that a change to lower water use land use in the GNN recharge zone can potentially increase recharge; 68% increase in the TLA and 16% to the TCSA. The effect of preferential flow on recharge to the unconfined aquifer and vertical leakage to the confined aquifers was unknown in the GNN recharge zone. In the absence of any other suitable techniques, we show that first order assessment of the recharge and vertical leakage can be assessed by sub-regional water balance estimates using a calibrated groundwater model. Management of the entire GNN recharge zone based on these new results is the key to ensuring the sustainability of the TCSA.

Acknowledgements

The editor and reviewers are thanked for evaluating the manuscript. We thank Glyn Ashman for manuscript review.

NOTES

^*Corresponding author.