Parameterization of a National Groundwater Model for New Zealand

: Groundwater is a vital source of water for humanity, with up to 50% of global drinking water and 43% of irrigation water being derived from such sources. Quantitative assessment and accounting of groundwater is essential to ensure its sustainable management and use. TopNet-GW is a parsimonious groundwater model that was developed to provide groundwater simulation at national, regional, and local scales across New Zealand. At a national scale, the model can help local government authorities estimate groundwater resource reliability within and between regions. However, as many catchments are ungauged, the model cannot be calibrated to local conditions against observed data. This paper, therefore, describes a method to derive an a priori, reach-scale groundwater model parameter set from national-scale hydrogeological datasets. The parameter set includes coefﬁcients of lateral (k l ) and vertical (k r ) conductivity and effective aquifer storage (S). When the parameter set was used with the TopNet-GW model in the Wairau catchment in the Marlborough region (South Island, New Zealand), it produced a poor representation of peak river ﬂows but a more accurate representation of low ﬂows (overall NSE 0.64). The model performance decreased in the smaller Opawa catchment (NSE 0.39). It is concluded that the developed a priori parameter set can be used to provide national groundwater modeling capability in ungauged catchments but should be used with caution, and model performance would beneﬁt greatly from local scale calibration. The parameter derivation method is repeatable globally if analogous hydrological and geological information is available and thus provides a basis for the parameterization of groundwater models in ungauged catchments. Future research will assess the spatial variability of parameter performance at a national scale in New Zealand.


Introduction
Accurate knowledge and accounting of groundwater resources are critical for effective water resource management, particularly in groundwater-dominated catchments.A range of groundwater model types are available, but parameterization of such models can be challenging where there are limited observed data available.Whilst the adoption of parsimonious models in some way eases this problem, a logical and transparent strategy for model parameter identification is still required [1][2][3].
The TopNet-GW model [4] is a parsimonious groundwater model based on the TopNet model, which is a national-scale rainfall-runoff model for New Zealand [5,6].The TopNet model performs best in large, medium-wet catchments but is insensitive to groundwater contribution [7,8].TopNet-GW was, therefore, developed to represent the groundwater contribution to surface water [4,9].Both TopNet and TopNet-GW models form part of the New Zealand Water Model (NZWaM) modeling framework [10].
By providing a greater representation of groundwater dynamics and the interaction of groundwater with surface water, TopNet-GW can help government authorities meet the requirements of New Zealand's National Policy Statement for Freshwater Management (NPS-FM), which requires quantitative modeling of groundwater resources.
Whilst commercial and research-driven models of aquifer systems have previously been created for New Zealand [11], different scales of application and underlying assumptions make direct comparison of results difficult.For example, recent models have included the use of eigenvectors [12], model ensemble smoothing techniques [13], water age distributions [14], three-dimensional modeling [15], and parameter up-scaling techniques [16].In addition, the use of multiple data aggregation methods on hydrological, geological, and soil input data (such as the River Environment Classification (REC) database, QMAP, and S-map [17][18][19]) can create a large range of structural variations in each of model type.

Aims of Study
Current adoption of the TopNet-GW model is hampered by a lack of knowledge about which parameter values should be used in the model across New Zealand.A method to define land cover, soil, and hydrological parameters from national datasets was developed by Bandaragoda et al. during the development of the TopNet surface water model [5], but no allowance was made for groundwater parameters.The primary aim of this study, therefore, is to develop a nationally consistent dataset of groundwater parameters, both for the TopNet-GW model and for potential use in other groundwater models.

Materials and Methods
The TopNet-GW model uses tight coupling surface and groundwater processes, which it achieves by passing water between surface and sub-surface stores within a single timestep.Whilst tightly coupled surface-groundwater models have been defined previously [20][21][22][23][24][25], none have been applicable at a national scale or used a priori parameterization.
TopNet-GW calculates lateral and vertical groundwater flux from groundwater storage capacity, hydraulic conductivity, river flow, and river channel dimensions.Flow (q) from any sub-catchment consists of three components: lateral flow to up to two downstream or adjacent groundwater sub-catchments (q l ); vertical flow to (or from) the river channel (q r ); and vertical flow to deeper groundwater (q x )): q = q l + q r + q x (1) Figure 1 illustrates the direction of flows in three adjacent groundwater sub-catchment stores (S 1 , S 2 , and S 3 ), where q l is the lateral flow between groundwater sub-catchments, q r is the flow between shallow groundwater and the river, and q x is the vertical loss from the groundwater sub-catchment store.
of groundwater with surface water, TopNet-GW can help government authorities mee the requirements of New Zealand's National Policy Statement for Freshwater Manage ment (NPS-FM), which requires quantitative modeling of groundwater resources.
Whilst commercial and research-driven models of aquifer systems have previously been created for New Zealand [11], different scales of application and underlying assump tions make direct comparison of results difficult.For example, recent models have in cluded the use of eigenvectors [12], model ensemble smoothing techniques [13], water age distributions [14], three-dimensional modeling [15], and parameter up-scaling technique [16].In addition, the use of multiple data aggregation methods on hydrological, geologi cal, and soil input data (such as the River Environment Classification (REC) database QMAP, and S-map [17][18][19]) can create a large range of structural variations in each o model type.

Aims of Study
Current adoption of the TopNet-GW model is hampered by a lack of knowledge about which parameter values should be used in the model across New Zealand.A method to define land cover, soil, and hydrological parameters from national datasets wa developed by Bandaragoda et al. during the development of the TopNet surface wate model [5], but no allowance was made for groundwater parameters.The primary aim o this study, therefore, is to develop a nationally consistent dataset of groundwater param eters, both for the TopNet-GW model and for potential use in other groundwater models

Materials and Methods
The TopNet-GW model uses tight coupling surface and groundwater processes which it achieves by passing water between surface and sub-surface stores within a single timestep.Whilst tightly coupled surface-groundwater models have been defined previ ously [20][21][22][23][24][25], none have been applicable at a national scale or used a priori parameteriza tion.
TopNet-GW calculates lateral and vertical groundwater flux from groundwater stor age capacity, hydraulic conductivity, river flow, and river channel dimensions.Flow (q from any sub-catchment consists of three components: lateral flow to up to two down stream or adjacent groundwater sub-catchments (ql); vertical flow to (or from) the rive channel (qr); and vertical flow to deeper groundwater (qx)): q = q l + q r + q x (1) Figure 1 illustrates the direction of flows in three adjacent groundwater sub-catch ment stores (S1, S2, and S3), where ql is the lateral flow between groundwater sub-catch ments, qr is the flow between shallow groundwater and the river, and qx is the vertical los from the groundwater sub-catchment store.Each outflow component from a groundwater sub-catchment (q l , q x , and q r , each with units of m/s) is assumed to be a linear function of groundwater stored in the sub-catchment (S, with units of m) such that: where k L is a lateral discharge coefficient [s where k R is the river discharge coefficient [s −1 ], and where k X is the deep groundwater percolation coefficient [s −1 ].As a result, outflow from each groundwater sub-catchment (q) can be described by a one-dimensional flow vector and is a function of stored groundwater (S) such that: where K lrx is a groundwater discharge coefficient (s −1 ) of the sub-catchment.It follows that: and thus To simplify further, we can also describe each flow coefficient ratio as a single variable (f): and We introduce the flow coefficient ratios because, with the limited information available, it is more practical to estimate a single outflow coefficient (K) and some ratios of flows (f) rather than several outflow coefficients.

Parameterization Scheme
National-scale datasets used to parameterize the above model include the New Zealand Groundwater Atlas [26] and water-table map [27][28][29]; QMAP 1:250,000 national geological map of New Zealand [13,30]; the New Zealand Aquifer Map [31]; and the NZ River Maps hydrology, ecology, and water quality metrics [32].Derived parameters from these coverages include depth to groundwater, porosity, and hydraulic conductivity [11,26,27]; river-bed dimensions, sediment type, and river-bed hydraulic conductivity; and surface-groundwater interaction [11].Figure 2 illustrates steps in the parameterization schema and datasets used.river-bed dimensions, sediment type, and river-bed hydraulic conductivity; and surfacegroundwater interaction [11].Figure 2 illustrates steps in the parameterization schema and datasets used.The national water-table map [27][28][29] also describes static water-table head (from which flow direction can be derived), sub-catchment scale saturated hydraulic conductivity (Ksat), and groundwater storage (S).River-bed hydraulic conductivity and channel width [32] were used to determine the river discharge coefficient (kr) and lateral flow coefficient (kl) for each reach.
Estimation of parameter sets was performed for each regional administration boundary to allow a national-scale model to be run by the region and to promote model adoption by local government agencies.

Surface-Groundwater Interaction (qr)
There are three potential methods to estimate flow between the surface and groundwater sub-catchments: Flow rate can be described as the average fraction of river flow to groundwater, such that: The national water-table map [27][28][29] also describes static water-table head (from which flow direction can be derived), sub-catchment scale saturated hydraulic conductivity (K sat ), and groundwater storage (S).River-bed hydraulic conductivity and channel width [32] were used to determine the river discharge coefficient (k r ) and lateral flow coefficient (k l ) for each reach.
Estimation of parameter sets was performed for each regional administration boundary to allow a national-scale model to be run by the region and to promote model adoption by local government agencies.

Surface-Groundwater Interaction (q r )
There are three potential methods to estimate flow between the surface and groundwater sub-catchments: Flow rate can be described as the average fraction of river flow to groundwater, such that: q r = − f r q sf (14) where f l ∈ [0, 1], q sf is the surface water flow, and f r can be determined from long-term river flow accretion data (i.e., variation in q r relative to q sf ).
1. NZ River Maps [32] coverage can be used to say that: q r = −W L K bed (15) where W and L are river-bed width and length, and K bed is the hydraulic conductivity of the river-bed [12].

2.
Like the above method but uses river depth such that: q r = − g(river) W L K bed (16) where g(river) is the fraction of river depth lost to groundwater per timestep.

Lateral Groundwater Flow (q l )
The lateral groundwater flow coefficient (f l ) is simplified such that: f l = 1 when water flows from the river to groundwater ('losing reach'); f l = 0 when water flows from groundwater to the river ('gaining reach').
However, f l can be given a value between 0 and 1 if prior information is available; thus, f l ∈ [0, 1].

Flow to Deeper Groundwater (q x )
Groundwater loss (to a non-returnable deeper groundwater zone) can be estimated using a simple water balance equation.Knowledge of other hydrological fluxes (rainfall, runoff, and evapotranspiration) can be used to calculate this term if losses are expected and there is insufficient reliable hydrogeological information available.Under most circumstances, a 'closed catchment' is assumed (i.e., f x = 0) unless evidence suggests otherwise.

Flow Coefficients
The sum of the flow coefficients (Equation ( 6)) can be estimated from hydraulic conductivity and porosity data derived from the New Zealand groundwater Atlas [26,27].Aquifer near-surface hydraulic conductivity (K) is defined [29,33] and calculated using: where K 0 is the hydraulic conductivity at or near the surface, z is the depth, and c is a function of terrain slope, climate, geology derived from mechanical and chemical denudation, and tectonic uplift rates of large sedimentary basins defined by the following equation: where s is the terrain slope and a, b, and c min are constants set to 75, 150, and 4, respectively [30].The resulting K sat values range from 0 to 87 m/day.Porosity (ϕ z ) at depth (z) is derived using the relationship described for sandstone and shales with clay dependence [34] such that: where z is depth, and ϕ 0 is effective porosity at the surface.Cl is the percentage clay factor, resulting in porosity ranging from 0 to 32%.Both K sat and storage (S) (calculated from porosity and estimated aquifer volume) are assumed to decrease with depth (z) below the surface.If we assume that the aquifer behaves like a linear reservoir, the characteristic hydraulic response time will be 1/K, so that: If the equilibrium water table is at depth z, then hydraulic conductivity and porosity are evaluated with respect to that depth.As both K sat (z) and S(z) are available at 200 m grid resolution of the water-table map, K can be calculated for individual TopNet-GW sub-catchment modeling units (i.e., Strahler order 1 or 3 catchments of approximately 0.5 km 2 and approximately 10 km 2 , respectively).
If K is known, then only two out of the three flow ratio parameters f l , f r , and f x described in Equations ( 10)-( 12) need to be estimated, as the third parameter can be calculated from Equation (13).As it is generally easier to obtain estimates of flow between the river and groundwater store (for example, through streamflow accretion survey or borehole monitoring) compared to an estimation of flow to the deeper groundwater store, f l and f r are estimated rather than f x .
For losing reaches (that exhibit a nett transfer of water from the river to groundwater), flow from shallow groundwater to the river (k r ) is zero, but there may be a lateral flow between each shallow groundwater catchment so that: For gaining reaches (that exhibit a nett transfer of water from groundwater to river), there will be a higher flow per unit area than is suggested by the larger aquifer water balance as not all the area of the aquifer is producing flow into the river.If a l and a g are the areas of the catchment with losing and gaining reaches, respectively, then we can say that: k l = K q l q (24)

Regional Geology
Figures 3 and 4 illustrate the distribution of lithology type within aquifers in each administrative region of New Zealand.Similarities in aquifer lithology between regions [35] suggest that model parameterization can be simplified by using combined region parameter sets.Northland is characterized by sand, silt and gravel, and volcanic aquifers.Surfacegroundwater interaction increases around the perimeter of the basalt and greywacke formations.Groundwater recharge typically occurs in upland areas, which then drain via springs in the foothills and to the lower valleys.Aquifers in the region are characterized as either quaternary coastal aquifers or inland volcanic aquifers.Auckland, to the south, differs slightly in that most groundwater abstraction is from the geothermal and volcanic basalt aquifers.MODFLOW and MIKE-SHE have been used [36] to model the hydrogeological response of the unconfined fractured basalt and semi-confined sandstone and mudstone aquifer systems.
The Waikato region can be divided geologically into three formations (alluvial deposits, volcanic, and sedimentary).The location and extent of groundwater availability are reflected in the location of consented groundwater takes, which are predominantly found in alluvial deposits [37].The Bay of Plenty region, which is dominated by rhyolite and ignimbrite lithology, with gravel, sand, and silts aquifers towards the coast, can be divided into three major groundwater catchments (Rangitikei, Tarawera, and Whakatāne).Significant parts of this region have previously been modeled by [38][39][40][41].Other models include a MODFLOW transient model for the Kaituna water management area [42], a flow and transport model for Tarawera [43], and a Western Bay of Plenty groundwater flow model [44].

Tasman and Nelson Marlborough Canterbury
West Coast Otago Southland      The aquifers within the Gisborne region surround Poverty Bay Flats, where over 1400 bores draw water for irrigation and domestic use from the gravel, silt, and sand aquifers [45,46].Except for the Te Hapara aquifer, none of the known sand aquifers have a direct connection to the sea [46].It is worth noting that 51% of aquifers in the region are 'limestone', which feeds into gravel flats or flows directly to the coast.All rivers close to the coast are likely to recharge groundwater in alluvial valleys and sand dune systems.
Hawke's Bay region exhibits similar aquifer systems (Heretaunga Plains and Ruataniwha Basin).The former is an alluvial plain formed by sediments deposited by the Ngaruroro, Tukituki, and T ūtaekurī Rivers [47].The Waipawa River through the Ruataniwha Plains indicates losses to groundwater of around 2.2 m 3 /s [48][49][50][51][52].With 81% of the defined aquifers in the region classed as 'limestone', rivers generally lose water to the aquifer upstream and gain water downstream.
Taranaki is different from all other North Island regions as it is dominated by andesitic volcanic cones.It has five principal aquifer systems named after the geological formations (Matemateāonga, Whenuakura, Tāngāhoe, Marine Terraces, and Taranaki Volcanics).Shallow groundwater is extracted from volcanic deposits and marine terraces east of New Plymouth and along the SW coast towards Wanganui [53], whereas deeper groundwater (>200 m bgl) is extracted from Tertiary sediments in north Taranaki.
The Upper Hutt aquifer in the Wellington region receives recharge from the Hutt River between Maoribank Bend and 700 m upstream of the Whakatikei confluence [30].Data suggest that 500-600 L/s is lost between Maoribank Bend and Pine Ave [46].The river crosses the Wellington Fault at Moonshine Bridge, where the riparian geology changes from low-permeability Pleistocene fan and alluvium deposits to more-permeable Holocene alluvium.The river provides significant recharge to the unconfined alluvial gravel aquifer upstream of the Whakatikei River confluence [54][55][56][57][58].Other significant aquifer bodies in the region include the Wairarapa [59][60][61] and Kapiti Coast [62].
Like the Wellington region, the general direction of groundwater movement in the Horizons (Manawat ū-Whanganui) region is from the inland hills towards the coast.The strata of the region consist of three main sequences: i. the geological basement made up of low permeability and heavily inundated greywacke (Ruahine and Tararua Ranges); ii.fine-grained marine sedimentary strata comprised of low permeability siltstone and/or mudstone; and iii.alluvial deposits formed by the erosion of the greywacke ranges.Thirtyeight percent of the region's aquifers are classed as limestone, with the remainder being gravel, sand, and silts.
In the South Island, virtually all aquifers are classed as 'gravel, silt or sand'.Whilst this should make parameterization easier, the 'gravel, silt or sand' classification is large, so refinement of derived model parameters may be required at the sub-catchment scale.

Results
Figure 5 illustrates the similarity between neighboring regions of Northland and Auckland; Waikato and Bay of Plenty; Gisborne and Hawkes Bay; and Wellington and Manawat ū-Whanganui.The extent of these similarities and the case for combined region parameter sets is described below.
Additional quantitative analysis of known parameter ranges (e.g., K sat and porosity) can aid in the estimation of derived parameter sets.Figure 6a,b, for example, illustrate the distribution of near-surface hydraulic conductivity and porosity.
Standardized values of each parameter are plotted to illustrate the variability of each parameter relative to rock type and region.Correlation between the hydraulic conductivity and porosity were calculated for each rock type and region, with strongest correlations occurring in basalt (R 2 = 0.97), basalt/breccia/scoria (R 2 = 0.93), and pumice/mudstone/sandstone (R 2 = 0.87) rock types, and in Taranaki (R 2 = 0.87), Gisborne (R 2 = 0.63); Southland (R 2 = 0.59), and Bay of Plenty (R 2 = 0.59) regions.Notably, there is limited correlation in the South Island regions, which, despite having a single geology type (gravel, sand, silt), have very large variability in values of both hydraulic conductivity and porosity.This leads to greater uncertainty in parameterization in those regions.
From the above analysis and description, it can be seen that adjacent regions may show some similarity in derived model parameter values (K sat and porosity) (Figure 6c).This happens because the dominant geology of any region, which is a primary control of ground surface parameters, often runs across the administrative region boundaries.This can also be seen at a national scale, as the variability of surface parameters is defined by underlying aquifer geology (Figure 6d).Regions with the most similar geology, and hence parameters, include Northland-Auckland-Waikato; Bay of Plenty-Gisborne; Hawke's Bay-Taranaki-Manawat ū-Whanganui; Wellington-Marlborough; Tasman-West Coast-Canterbury; and Otago-Southland.

Results
Figure 5 illustrates the similarity between neighboring regions of Northland and Auckland; Waikato and Bay of Plenty; Gisborne and Hawkes Bay; and Wellington and Manawatū-Whanganui.The extent of these similarities and the case for combined region parameter sets is described below.
Additional quantitative analysis of known parameter ranges (e.g., Ksat and porosity) can aid in the estimation of derived parameter sets.Figure 6a Standardized values of each parameter are plotted to illustrate the variability of each parameter relative to rock type and region.Correlation between the hydraulic conductivity and porosity were calculated for each rock type and region, with strongest correlations occurring in basalt (R 2 = 0.97), basalt/breccia/scoria (R 2 = 0.93), and pumice/mudstone/sandstone (R 2 = 0.87) rock types, and in Taranaki (R 2 = 0.87), Gisborne (R 2 = 0.63); Southland (R 2 = 0.59), and Bay of Plenty (R 2 = 0.59) regions.Notably, there is limited correlation in the South Island regions, which, despite having a single geology type (gravel, sand, silt), have very large variability in values of both hydraulic conductivity and porosity.This leads to greater uncertainty in parameterization in those regions.
From the above analysis and description, it can be seen that adjacent regions may show some similarity in derived model parameter values (Ksat and porosity) (Figure 6c).This happens because the dominant geology of any region, which is a primary control of ground surface parameters, often runs across the administrative region boundaries.This can also be seen at a national scale, as the variability of surface parameters is defined by underlying aquifer geology (Figure 6d).Regions with the most similar geology, and hence parameters, include Northland-Auckland-Waikato; Bay of Plenty-Gisborne; Hawke's Bay-Taranaki-Manawatū-Whanganui; Wellington-Marlborough; Tasman-West Coast-Canterbury; and Otago-Southland.
The resulting vertical (kr) and lateral (kl) hydraulic conductivity coefficients for losing and gaining reaches (as defined in Equations ( 21)-( 24)) are illustrated in Figure 7.The distribution of the lateral flow coefficient (kl) more strongly reflects the variation in geological and land surface parameters, whereas the vertical coefficient is controlled by the extent of surface-groundwater interaction.The resulting vertical (k r ) and lateral (k l ) hydraulic conductivity coefficients for losing and gaining reaches (as defined in Equations ( 21)-( 24)) are illustrated in Figure 7.The distribution of the lateral flow coefficient (k l ) more strongly reflects the variation in geological and land surface parameters, whereas the vertical coefficient is controlled by the extent of surface-groundwater interaction.

Uncertainty
Whilst regional councils in New Zealand have both surface and groundwater monitoring networks, there remains much uncertainty in the total extent and magnitude of surface-groundwater interaction [63].Surface-groundwater interaction in New Zealand

Uncertainty
Whilst regional councils in New Zealand have both surface and groundwater monitoring networks, there remains much uncertainty in the total extent and magnitude of surface-groundwater interaction [63].Surface-groundwater interaction in New Zealand predominantly occurs after rivers have emerged from the high mountain ranges and as they flow across alluvial plains.River flow depletion and accretion surveys can be used to estimate the magnitude of surface-groundwater flow, but this tends to be spatially and temporally dynamic, as it is dependent on local geological variability and sub-surface hydraulic connectivity [64].
The parameter sets outlined in this study inherit assumptions contained within the hydrogeological datasets from which they are derived.For example, the surface-groundwater interaction map [9] is based on the assumption of three characteristic landscape types, including permeable headwaters (accreting rivers), inland flood plains (depleting rivers), and accreting rivers (coastal catchments).As a result, only sub-catchments within these landscape types are defined as active groundwater catchments within the parameter datasets.
Simulated river flow from the TopNet-GW model has been shown to be most sensitive to hydraulic conductivity between the river and aquifer (k r ) [9].However, simulations are also sensitive to the extent of shallow aquifer zone and lateral groundwater flow (k l ).Thus, whilst we would expect river flow to be most sensitive to k r , it will be mitigated by k l to an extent dependent on the existence of a shallow aquifer layer and the aquifer storage potential (represented by porosity).
The uncertainty associated with the parameter estimation method used is related to the water-table depth map [30], which assumes a state of equilibrium within modeled aquifers.The implication of this is that the derived flow coefficients k r and k l underestimate flow gradients during model simulations of rapid groundwater depletion or recharge.The impact of this uncertainty can be reduced through local-scale calibration and validation through field study.Similarly, the derivation of depth to hydraulic basement value was based on the analysis of large-scale fluvial sedimentary basins [30].Thus, the estimates of sub-surface hydraulic conductivity made for other landform types (e.g., small coastal basins or volcanic aquifers) are susceptible to greater uncertainty (as estimates were based on hydraulic depth to basin).The authors also noted that no consideration of the impact of the potential of confining layers or material anisotropy on hydraulic conductivity was made.

Case Study: Marlborough, NZ
The above datasets were used to parameterize TopNet-GW for the Marlborough region (South Island).The locations of river catchments that exhibit surface-groundwater interaction [9] are shown in Figure 8a,b.The mean annual groundwater head (m above mean sea level) and depth to water table (m below ground level) (as defined by the watertable map [27]) are shown in Figure 8c,d.Dominant groundwater flow direction is derived from the distribution of mean annual groundwater head (Figure 8e), and the location of gauged streams in Marlborough is shown in Figure 8f.
Figure 9a,b illustrate distributions of conductivity (K) and aquifer storage (S) for the region, as derived from the water-table height, river-bed conductivity, and channel width (Figure 9c,d).
Figure 9e,f show coefficients of vertical and lateral flow (k r and k l ).It can be seen that k r and river-bed conductivity are spatially correlated as river-bed permeability is a direct control of surface-groundwater interaction through the river channel.
To assess the utility of the above parameter sets for groundwater modeling, TopNet-GW was run for the Marlborough region for the period 2000 to 2010.The performance of the TopNet_GW model was then compared to the TopNet surface water model and observed data in the Wairau River and Opawa River (see Figure 10).Both models performed well in the larger Wairau River catchment (Nash-Suttcliff Efficiency (NSE)-0.64).Performance was worse in the smaller Opawa River (NSE = 0.3).
The flow duration curves and cumulative flow diagrams in Figure 10 illustrate how each model replicates annually averaged flow characteristics.Both TopNet-GW and TopNet models underestimate lower quartile flows in the Wairau River.This indicates that the models are over-predicting 'leakage' from the river to the groundwater store.By contrast, in the Opawa River catchment (modeled from mid-2004), both models overestimate peak flows.The TopNet-GW model performs better than the TopNet model, which underestimates low flows to a greater extent.

Case Study: Marlborough, NZ
The above datasets were used to parameterize TopNet-GW for the Marlborough region (South Island).The locations of river catchments that exhibit surface-groundwater interaction [9] are shown in Figure 8a,b.The mean annual groundwater head (m above mean sea level) and depth to water table (m below ground level) (as defined by the water-table map [27]) are shown in Figure 8c,d.Dominant groundwater flow direction is derived from the distribution of mean annual groundwater head (Figure 8e), and the location of gauged streams in Marlborough is shown in Figure 8f.
Figure 9a,b illustrate distributions of conductivity (K) and aquifer storage (S) for the region, as derived from the water-table height, river-bed conductivity, and channel width (Figure 9c,d).x FOR PEER REVIEW 14 of 20 (e) (f) Figure 9e,f show coefficients of vertical and lateral flow (kr and kl).It can be seen that kr and river-bed conductivity are spatially correlated as river-bed permeability is a direct control of surface-groundwater interaction through the river channel.
To assess the utility of the above parameter sets for groundwater modeling, TopNet-GW was run for the Marlborough region for the period 2000 to 2010.The performance of the TopNet_GW model was then compared to the TopNet surface water model and observed data in the Wairau River and Opawa River (see Figure 10).Both models performed well in the larger Wairau River catchment (Nash-Suttcliff Efficiency (NSE)-0.64).Performance was worse in the smaller Opawa River (NSE = 0.3).
The flow duration curves and cumulative flow diagrams in Figure 10 illustrate how each model replicates annually averaged flow characteristics.Both TopNet-GW and

Discussion
The parameterization scheme described above is used to produce a reach-scale groundwater-model parameter set for aquifers in New Zealand.Vertical and lateral hydraulic conductivity and aquifer storage are defined from hydrogeological data.Whilst the method is globally reproducible, there are several features in the approach that are dependent on the availability of key datasets.
For example, the method uses a high-resolution digital elevation model (DEM) to define a hydrologically continuous drainage network.National coverages of geological data, aquifer locations, and surface water catchment boundaries, which are available in most countries, are then defined relative to the DEM.By contrast, river channel characteristics, river-bed conductivity, and estimates of depth to groundwater may be more difficult to identify at the required resolution.
Comparison of model parameters and output relative to different aquifers in New Zealand is likely to yield varying results, as local or regional scale aquifer models may be of varying resolution, extent, and spatial variability.Whilst outside the scope of this study, it is envisaged that such a comparison will help inform the refinement of model parameter sets in the future.
In the case study for the Marlborough region, the performance of the parameter set, as represented by the TopNet-GW model output, is regionally specific, and results cannot be assumed to be representative of the whole country.Future research will assess the spatial variability of model output along river profiles within aquifer areas and the seasonal variability in surface-groundwater interaction, as determined by water-table fluctuation, and compare model internal variables with observed groundwater height variation.

Discussion
The parameterization scheme described above is used to produce a reach-scale groundwater-model parameter set for aquifers in New Zealand.Vertical and lateral hydraulic conductivity and aquifer storage are defined from hydrogeological data.Whilst the method is globally reproducible, there are several features in the approach that are dependent on the availability of key datasets.
For example, the method uses a high-resolution digital elevation model (DEM) to define a hydrologically continuous drainage network.National coverages of geological data, aquifer locations, and surface water catchment boundaries, which are available in most countries, are then defined relative to the DEM.By contrast, river channel characteristics, river-bed conductivity, and estimates of depth to groundwater may be more difficult to identify at the required resolution.
Comparison of model parameters and output relative to different aquifers in New Zealand is likely to yield varying results, as local or regional scale aquifer models may be of varying resolution, extent, and spatial variability.Whilst outside the scope of this study, it is envisaged that such a comparison will help inform the refinement of model parameter sets in the future.
In the case study for the Marlborough region, the performance of the parameter set, as represented by the TopNet-GW model output, is regionally specific, and results cannot be assumed to be representative of the whole country.Future research will assess the spatial variability of model output along river profiles within aquifer areas and the seasonal variability in surface-groundwater interaction, as determined by water-table fluctuation, and compare model internal variables with observed groundwater height variation.
A significant drawback of the existing parameter set is that it has been developed with reference to a steady-state surface-groundwater interaction map.As a result, modeled groundwater flow is unidirectional, either as flow from the river to the groundwater store or as flow from the groundwater store to the river.This reduction in model complexity results in lower hydrograph peaks and elevated low-flow characteristics.
The results of this research build on the conceptual understanding of previous work that identified potential groundwater recharge zones.Whilst there is some overlap in methodology with the previous study (same land-use, soil, geology, and topographic data), both studies identify the potential and relative magnitude of groundwater recharge.In this study, the component datasets were further aggregated to provide a reach-scale characterization of subsurface properties.

Conclusions
An a priori parameterization method was developed for the TopNet-GW model that allowed it to be applied in ungauged catchments in New Zealand.A nationwide parameter set was produced to represent regional aquifer storage behavior and groundwater flux.When the model was applied in the Marlborough region, it performed well in predicting low flows but less well in predicting peak flows.The performance of the model indicated that the spatial distribution of groundwater processes is correct but that local scale calibration may be required to improve model accuracy.
Improved a priori parameter identification may be achieved by using additional hydrogeological data (such as the aquitard and aquiclude layers and their associated properties).However, in the first instance, a more thorough assessment of the spatial variation in model performance and parameter uncertainty across all regions of New Zealand is required.

Figure 1 .
Figure 1.Conceptual diagram of flows from the groundwater stores (sub-catchments) of th TopNet-GW model.

Figure 1 .
Figure 1.Conceptual diagram of flows from the groundwater stores (sub-catchments) of the TopNet-GW model.

Figure 2 .
Figure 2. Schematic diagram relating database sources (right) to parameters required in groundwater modeling.

Figure 2 .
Figure 2. Schematic diagram relating database sources (right) to parameters required in groundwater modeling.

Figure 3 .
Figure 3. North Island aquifer geology by region (not to scale); color key is shown in Figure 5.

Figure 3 .
Figure 3. North Island aquifer geology by region (not to scale); color key is shown in Figure 5.

Figure 4 .
Figure 4. South Island aquifer geology by region (not to scale); color key is shown in Figure 5.

Figure 5 .
Figure 5. Percentage area of aquifer rock type within each region (red boxes circling similar regions).

Figure 4 .
Figure 4. South Island aquifer geology by region (not to scale); color key is shown in Figure 5.

Figure 4 .
Figure 4. South Island aquifer geology by region (not to scale); color key is shown in Figure 5.

Figure 5 .
Figure 5. Percentage area of aquifer rock type within each region (red boxes circling similar regions).Figure 5. Percentage area of aquifer rock type within each region (red boxes circling similar regions).

Figure 5 .
Figure 5. Percentage area of aquifer rock type within each region (red boxes circling similar regions).Figure 5. Percentage area of aquifer rock type within each region (red boxes circling similar regions).

Figure 6 .Figure 6 .
Figure5illustrates the similarity between neighboring regions of Northland and Auckland; Waikato and Bay of Plenty; Gisborne and Hawkes Bay; and Wellington and Manawatū-Whanganui.The extent of these similarities and the case for combined region parameter sets is described below.Additional quantitative analysis of known parameter ranges (e.g., Ksat and porosity) can aid in the estimation of derived parameter sets.Figure6a,b, for example, illustrate the distribution of near-surface hydraulic conductivity and porosity.

Figure 6 .
Figure 6.Spatial distribution of aquifer near-surface: (a) K sat (m/day) and (b) porosity (%); boxplots of standardized median, inter-quartile range, maxima, minima, and single outlier values by (c) region (with red boxes indicating similar regions); and (d) rock type.

Figure 7 .
Figure 7. Spatial distribution of (a) derived kl parameters and (b) derived kr parameters.

Figure 7 .
Figure 7. Spatial distribution of (a) derived k l parameters and (b) derived k r parameters.

Figure 10 .
Figure 10.Cumulative flow and flow duration curve of simulated and observed streamflow on the Wairau (Barnetts Bank) and Opawa (at Blicks Lane).

Figure 10 .
Figure 10.Cumulative flow and flow duration curve of simulated and observed streamflow on the Wairau (Barnetts Bank) and Opawa (at Blicks Lane).