Application of global models and satellite data for smaller-scale groundwater recharge studies

Large-scale models and satellite data are increasingly used to characterise groundwater and its recharge at the global scale. Although these models have the potential to fill in data gaps and solve trans-boundary issues, they are often neglected in smaller-scale studies, since data are often coarse or uncertain. Large-scale models and satellite data could play a more important role in smaller-scale (i.e., national or regional) studies, if they could be adjusted to fit that scale. In New 5 Zealand, large-scale models and satellite data are not used for groundwater recharge estimation at the national scale, since regional councils (i.e., the water managers) have varying water policy and models are calibrated at the local scale. Also, some regions have many localised ground observations (but poor record coverage), whereas others are data-sparse. Therefore, estimation of recharge is inconsistent at the national scale. 10 This paper presents an approach to apply large-scale, global, models and satellite data to estimate rainfall recharge at the national to regional scale across New Zealand. We present a model, NGRM, that is largely inspired by the global-scale WaterGAP recharge model, but is improved and adjusted using national data. The NGRM model uses MODIS-derived ET and vegetation satellite data, and the available nation-wide datasets on rainfall, elevation, soil and geology. A valuable addition to the 15 recharge estimation is the model uncertainty estimate, based on variance, covariance and sensitivity of all input data components in the model environment. This research shows that, with minor model adjustments and use of improved input data, largescale models and satellite data can be used to derive rainfall recharge estimates, including their uncertainty, at the smaller scale, i.e., national and regional scale of New Zealand. The estimated New 20 Zealand recharge of the NGRM model compare well to most local and regional lysimeter data and recharge models. The NGRM is therefore assumed to be capable to fill in gaps in data-sparse areas and to create more consistency between datasets from different regions, i.e., to solve trans-boundary issues. This research also shows that smaller-scale recharge studies in New Zealand should include larger boundaries than only a (sub-)aquifer, and preferably the whole catchment. This research points 25 out the need for improved collaboration on the international to national to regional levels to further merge large-scale (global) models to smaller (i.e., national or regional) scales. Future research topics should, collaboratively, focus on: improvement of rainfall-runoff and snowmelt methods; inclusion of river recharge; further improvement of input data (rainfall, evapotranspiration, soil and geology); and the impact of recharge uncertainty in mountainous and irrigated areas. 30


Introduction
Large-scale hydrological models are increasingly used to characterise groundwater at the global scale (e.g., Gleeson et al., 2011;Fan et al., 2013), its fluxes (e.g., Döll and Fiedler, 2008;Beck et al., 2013) and its depletion (Wada et al., 2010).Furthermore, satellite data have been shown to be important for global-scale assessments of important water cycle variables, such as rainfall (GPM: 35 especially near or beyond the boundaries of study areas, which are often only a sub-region and not a whole catchment.Application of multiple local models in one area also shows inconsistency.For example, the impact of rainfall recharge uncertainty was demonstrated by White et al. (2003), who 75 used three rainfall recharge models for the Central Plains of Canterbury, New Zealand.They showed that estimated groundwater use as a percentage of rainfall recharge was highly model-dependent, and because of that dependency varied in between 63% and 80% in a relatively dry year.However, uncertainty has not typically been assessed by sub-regional models in New Zealand.Uncertainties can either be caused by the different model equations or by difference in input data.Most recharge 80 models are some form of soil water balance Rushton et al. (2006); Westenbroek et al. (2010); White et al. (2003), and only show difference in the derivation of actual ET, or the assumed rainfall-runoff ratios.A regional study in New Zealand by Rawlinson et al. (2015) showed that the difference in three recharge models was most likely caused by the inconsistency of model input data and not by the differences in model equations.Development of a large-scale model that estimates model 85 uncertainty and uses the same input data for the nation, is thus beneficial, even at the local scale.
This study describes the use of a national groundwater recharge model that is inspired by the global-scale WaterGap recharge model of Döll and Fiedler (2008) and uses MODIS satellite data of ET (Mu et al., 2011) and vegetation (Samanta et al., 2011).The aim of this research is to show that, with minor model adjustments and use of improved input data, large-scale models and satellite data 90 can be used to derive rainfall recharge estimates, including their uncertainty, on the national scale of New Zealand, or smaller.

Method
The model in this study is called the 'National Groundwater Recharge Model' (NGRM).We present the NGRM in three stages.First, a general description of the method is given.Second, input data are 95 elucidated, including a description of pre-processing steps.Third, detailed methodology, including uncertainty estimation, is described.

General description
The NGRM has a 1 km model grid that covers most of New Zealand (Figure 1).The model calculates rainfall recharge to groundwater in monthly time steps, currently covering the period of January 2000 100 to December 2014.The NGRM uses a simple soil water balance, i.e., it calculates a mass balance between vertical water inflow and outflow in a simple representation of the soil (i.e., one layer), and it is assumed that any surplus water in the soil layer either drains to groundwater or goes to runoff.
Other common soil water balance models are the groundwater recharge module in WaterGAP (Döll and Fiedler, 2008), USGS-SWB (Westenbroek et al., 2010), Rushton (Rushton et al., 2006).The 105 NGRM consists of a newly developed method and code, but is largely inspired by the approach of the WaterGAP model, that calculates an initial recharge and then uses several correction factors to estimate true recharge, i.e., for slope, soil, geology and more.The strengths of the WaterGAP model are, in our opinion: it is a simple but rigorous method; it is a well described and accepted method; and it includes a deeper 'geology' layer.The input data used for the NGRM model include 110 satellite data of ET and vegetation, and New Zealand national datasets of soil, geology, elevation, and precipitation.
The NGRM is uncalibrated.However, case study comparisons described within this article have been used to improve the nation-wide model equations such that the model is considered to provide suitable estimates in these case studies.

Input data and processing
The rainfall recharge model uses input data on precipitation, evapotranspiration, vegetation, topography, soil parameters, and geology.These are summarised below.All data were re-gridded to the monthly 1km x 1km model cells, which are projected to the national New Zealand Transverse Mercator (NZTM) grid.This was done by averaging (if the time step or grid resolution of the original 120 input data was smaller) or by its nearest value (if the grid resolution of the original input data was larger).

Precipitation
Daily precipitation data from 1972 to the present day are available for New Zealand in a regular grid (0.05 • of latitude and longitude, or approximately 5 km) from the Virtual Climate Station (VCS) 125 network (Tait, 2014).These data are described by Tait et al. (2006).Only data covering the period 1-Jan-2000 to 31-Dec-2013 were used in this study (mean annual values are shown in Figure 1a).

Evapotranspiration
The MOD16 algorithm uses satellite data from NASA's MODerate resolution Imaging Spectroradiometer (MODIS) sensor.The satellite-derived parameters embedded in the MOD16 dataset are 130 land cover, albedo, leaf area index (LAI), fraction of absorbed photosynthetically active radiation (FPAR), and enhanced vegetation index (EVI).Temperature, radiation, air humidity and pressure data are derived from the daily global meteorological reanalysis data set from NASA's Global Modeling and Assimilation Office (GMAO).The MOD16 algorithm, described in Mu et al. (2011), uses the Penman-Monteith approach to calculate evapotranspiration.The spatial resolution of the 135 MOD16 cells is approximately 1 km x 1 km.MOD16 PET and MOD16 AET data are available online in HDF (Hierarchical Data Format) files in 8-daily, monthly, and yearly intervals (NTSG, 2013).
Monthly MOD16 PET data were tested in the New Zealand setting by Westerhoff (2015), who also derived a MOD16 Penman PET conversion and uncertainty estimate for these data (since Penman Hydrol. Earth Syst. Sci. Discuss., doi:10.5194/hess-2016-410, 2016 Manuscript under review for journal Hydrol.Earth Syst.Sci.Published: 30 August 2016 c Author(s) 2016.CC-BY 3.0 License.
PET is the national standard).Westerhoff (2015) furthermore suggests that MOD16 AET could be

Terrain model
The Geographx New Zealand DEM 2.1 (Geographx, 2012) is a national digital elevation model (Figure 1c) on an 8 m x 8 m grid that was derived from a combination of New Zealand-based topographic data from Land Information New Zealand and data from the satellite-derived Shuttle Radar Topography Mission (SRTM: USGS, 2006).150

Geology
Geology is defined as the subsurface underlying the soil.The geology of the shallow subsurface, i.e. up to approximately 10 m depth, is described by the 1:250,000 geological map of New Zealand (QMAP: GNS Science, 2012).This digital map, a GIS polygon file, includes several polygon attributes, of which main rock type, secondary rock type, and age were used in this study.The main 155 rock type and secondary rock type were interpreted to hydrolithological classes and intrinsic permeability κ [m 2 ].The estimation of κ was based on a classification method described by Gleeson et al. (2011), who derive κ and standard deviation for 7 main hydrolithological classes, using the information from Freeze and Cherry (1979), their table 2.2, and calibration results of a multitude of global and local hydrological models.Additionally, this classification was put in the New Zealand context 160 by with 3 additional hydrolithological classes.These classes, each with a median and standard deviation value of κ (Table 1), are used in this research.All QMAP different rock type attributes could be summarised to 183 'dictionary' values of rock type, which were then interpreted to permeabilities with a look-up table approach.Two examples are given: (1) sandstone, greenstone were interpreted as 'coarse-grained sedimentary' and were given a permeability (log κ) of -12.5, with an uncertainty 165 (standard deviation 1σ) of 0.9; (2) limestone and shell beds were interpreted as carbonate and were given a permeability of -11.8, with an uncertainty of 1.5.
Permeabilities of main rock type and secondary rock type were combined by a weighting average, in which the κ value of main rock type was weighted two times heavier than the secondary rock type.

170
It was assumed that vertical K is equal to horizontal saturated hydraulic conductivity.For the conversion of intrinsic permeability to hydraulic conductivity K [m/day] we used Eq.1: following (Freeze and Cherry, 1979, their Eq. 2.28), where: µ is the dynamic viscosity of freshwater at 13 • C (=1.2155×10 3 kg /m s);
The estimation of decrease of K over age is not straightforward.For example, Akbar et al. (1995) point out, for carbonates, that the relation between porosity and age is influenced by many factors such as diagenesis, pressure, temperature, and erosion.We chose to include a correction factor with 180 age to exclude old (hard)rock, to show a more realistic distribution of aquifers throughout New Zealand.The correction factor is based on an exponential decrease of hydraulic conductivity with age (Eq.2): where T is the geological age [Ma] and α a constant that controls the decrease rate.The value of α 185 was (arbitrarily) set to 40, so that f T was close to 1 for Quaternary rock types (i.e., younger than 1.8 Ma).
A national map of estimated hydraulic conductivity as used for the NGRM is shown in Figure 1d.

Soil
Soil permeability is the rate that water moves through saturated soil.New Zealand-wide soil per-  2).The FSL was chosen because it is the only dataset that covers the entire nation.Soil permeability was interpreted from the original soil qualitative classes to relative soil permeability ratio f soil in between 0 195 and 1 (Table 2).Soil permeability ratios are shown in Figure 1e.
Profile of Available Water (PAW) is the difference between field capacity and permanent wilting point.PAW is a soil property defining the total amount of available water that can be held by a soil profile (Smith and Mullins, 2001).PAW data are available from the LRIS portal as an ArcGIS shapefile, with each soil class having a minimum and maximum value (Table 3).Modal values were 200 calculated according to Newsome et al. (2008) and shown in Figure 1f.For areas where PAW was unknown, the mean of all PAW modal values was assumed.

Leaf Area Index
Leaf Area Index (LAI) is derived from surface reflectance in the red and near-infrared bands measured by NASA's MODIS sensor and is defined as the one-sided green leaf area per unit ground area so that a conservative estimate of 2% of rainfall will be intercepted by the canopy at maximum LAI (approximately 6).The interception was subtracted from the precipitation dataset.It was assumed that the satellite-measured MOD16 AET already embeds the evaporation of the intercepted water and therefore Ic was not added to AET.

Snow
The monthly precipitation data, after correction for interception, were corrected for the precipitation being snow rather than rain.A very simple assumption was made: if the average monthly temperature (from the VCS data) is lower than 0, all precipitation is assumed snow.If the average monthly temperature is between 0 and 4, 30% of the precipitation is assumed snow.The implications of this 225 simplicity are further explained in the discussion.

Rainfall recharge model equation
In monthly time steps, rainfall recharge RRECH (in mm) was calculated for time step i as: where: 230 R i is surface rainfall (in mm, after correction for interception and snow); f slope is a correction factor for rainfall runoff due to slope [0 to 1]; AET is actual evapotranspiration (in mm); S is soil storage (mm); If RRECH i < 0 (Eq.4), then RRECH i is set to zero and the resulting storage deficit S i < S i−1 is used for the calculation of RRECH i+1 .
The slope correction factor f slope was calculated from the terrain model.After calculating a slope α (in degrees) from the terrain model, this was then used to calculate an initial runoff of the rainfall 240 due to slope: where the error function (erf) fits the empirical slope -runoff relation used by Döll and Fiedler (2008).The NGRM model prefers to use the simple relation to terrain slope, instead of other alternatives that use a relation between rainfall intensity and soil type (e.g., Rushton et al., 2006; 245 Westenbroek et al., 2010).Also, the NGRM model prefers the f slope to directly affect rainfall, and not affect AET and S, as in Döll and Fiedler (2008).Possible model limitations arising because of this simplification are explained in the discussion.
The initial soil storage S init for the first time step (January 2000) was estimated as a function of the modal values of PAW for a vegetated surface: This initial soil storage stems from the 'water holding capacity' of a soil with vegetation (i.e., roots) and is described by White et al. (2003) and Scott (2004) for a grass surface.Using this initial storage means the model has an unknown error in soil storage in the first few months of the model runs (because not all soils were at their driest in January 2000).It is assumed that this unknown error 255 is resolved after six months (July 2000), when most soils are at their wettest.
The factor f soil represents the impedance to recharge from soil and is given by the soil permeability ratio (Table 2).For f soil , it was assumed that soils lower than a value of 0.15 (equivalent to a soil permeability of less than 4mm/hr, Newsome et al., 2008) reject 75% of the rainfall recharge to surface runoff; all other soils accept all recharge.The model assumes an isotropic permeability, i.e., 260 that vertical and lateral flow are equal, and that any rejected recharge will add to runoff within the same model cell.
The geology factor f geol represents the impedance to recharge from geology.Hydraulic conductivity was assumed to only play a limiting role if its value was less than the rainfall recharge.In practice, this means that only very impermeable geology (e.g., clay, basement rock, mudstone) would limit 265 rainfall recharge.The factor f geol was then calculated as the ratio (in between 0 and 1) of potential rainfall recharge and hydraulic conductivity.

Uncertainty of rainfall recharge
The uncertainty of the NGRM model was calculated by error propagation, i.e., propagating the variance and covariance of all input components, including their model sensitivity, through the rainfall 270 recharge equation (Eq.4), after Tellinghuisen (2001): where σ f 2 is the variance of a function f , which has n = 1 : N input components; V is the variance-covariance matrix of all input components; g is a vector of input component ∂f /∂n i ; and g T is the transpose of g.Considering Eq. 4, the components of g are (RRECH is denoted as Ψ to 275 spare length of the equations): Where the components f soil and f geol were multiplied to one factor f soilgeol beforehand.The components of V are: The dots (i.e., '..') indicate covariances between the parameters.Covariance was calculated with a subset of 5% of all (100,000 of approximately 267,000 pixels and 2 of 14 years) monthly time series of all input components (R, f slope , AET, S, f soil and f geol ).The uncertainty of the model input components were estimated, as values of the standard deviation (1σ), as follows: -The uncertainty of rainfall per monthly time step was assumed 5% in plain areas (with many -The uncertainty in the estimation of f slope was assumed to be 10%.This was assumed to be affected by general inaccuracy of terrain models (e.g., Westerhoff et al., 2013) and averaging of the terrain model to the model grid;

295
-The uncertainty in daily Penman PET in New Zealand can be 10 -40% (Westerhoff, 2015) and is a function of the PET value.Daily values of this uncertainty function were compiled to monthly and mean annual values.For this study, it was assumed that the uncertainty of AET decreases with the AET/PET ratio; -Uncertainty in storage was assumed to be a function of PAW.A Gaussian distribution was 300 assumed between minimum and maximum values of PAW, making the 1σ value 16% of the range; -The standard deviation of the f soil values was assumed to be 10%, except for the 'Slow' and 'Rapid' classes, where they are chosen as 5%; -Uncertainty in f geol was assumed to be affected by the uncertainty of hydraulic conductivity.

305
Uncertainty in K can be very high: even higher than the actual value of K (Gleeson et al., 2011;Tschritter et al., 2014).This study assumes that the maximum standard deviation in K is less than or equal to K itself.The uncertainty was further assumed to only play a role if the recharge was higher than K (with both recharge and K converted to match the monthly time step).

310
Rescaling the covariance matrix to the values of the uncertainty in the input components was done as follows: where: σ 2 xy,scaled is the scaled variance;

Model equation limitations
The NGRM is considered a simplified model that aims for: a national and inter-regional overview; relating differences in existing local models; and estimating rainfall recharge in unexplored territory.
The modelled rainfall recharge and its uncertainty estimates are therefore also considered simplified, although recharge estimates and their uncertainty fit well to observed differences with measurements 410 and other local models.Because the model has not been calibrated on a smaller scale, like local models, it has to use generic or simplified assumptions.These assumptions are discussed below.
The slope-runoff relation is based on a sparse dataset of empirical values by Döll and Fiedler (2008), who relate total recharge to a slope correction factor.Other relations were explored, but none better were found for use on the national scale.Probably the best-known alternative is a curve 415 number (CN) approach relating soil type and land cover to runoff (Cronshey, 1986).This method has been developed for gently sloping hills in the United States and might not be able to deal with steeper slopes, e.g., New Zealand mountainous terrain.Other rainfall recharge models use a relation between rainfall intensity and soil (Rushton et al., 2006;Westenbroek et al., 2010) without taking into account terrain slope.The NGRM prefers to use the simple relation to terrain slope, because of 420 three reasons.First of all, rainfall intensity is not captured well on a monthly scale.Second, we are of the opinion that steeper slopes lead to rainfall 'splash' (Hendriks, 2010) and can form preferential flow channels, ultimately adding to surface runoff.Third, terrain slope is inherently related to soil.
For example, regardless of mean annual precipitation, long-term mean denudation rate in river basins is proportional to basin relief Ahnert (1970); Summerfield and Hulton (1994), leaving less erodible 425 material (e.g., soil) on the steep slopes, which amplifies the splash effect.
A simple snow correction was assumed to estimate rainfall from precipitation.This correction factor was based on a coarse assumption that precipitation is snow when temperature is below 0 • C.
A better defined snow-rainfall relation should be implemented in future updates of the NGRM.Another shortcoming of this simple correction is snowmelt: in Spring a substantial amount of the snow will melt, which is likely to substantially recharge the groundwater and have a large impact on the whole water budget White (2007).Further research should therefore be applied to implement a module of 'snow correction and snowmelt recharge'.
Heterogeneity and model up-scaling and down-scaling could cause a large, but unknown uncertainty, e.g., averaging high resolution soil and subsurface parameterisation over a 1 km grid cell.
Furthermore, averaging slope for 1 km grid cells could lead to a wrong estimate of runoff, which would loop back to wrong recharge and more uncertainty in the slope-runoff relation.A higher resolution representation of elevation and slope was not implemented at the national scale, since that would require significantly more computational power.For example, a typical run for 2000 to 2013 now takes 1 hour on a standard desktop computer.Using higher resolution data, such as 100 m, would make the input datasets 100 times larger, and parallelisation would need to be implemented for more efficient computation.However, smart solutions can be explored, e.g., by compiling the high resolution only as spatial statistics of a 1 km x 1 km pixel, instead of using the full high resolution data arrays.Future research on application of the NGRM model on the local scale should therefore address these issues of scaling and heterogeneity.
Hydraulic conductivities of the underlying geology in this research are given for saturated flow.K values for unsaturated flow can be lower than for saturated flow, but they are non-linear and not easy to estimate Hendriks (2010); Fitts (2013).The simple NGRM model does not calculate unsaturated K vales for several reasons.First, since the uncertainty in K is already high and was therefore in the model equation already clipped to the actual value of K. Second, it was considered to adjust K relative to soil water deficit by simply choosing a lower value of saturated K, e.g., 75% of saturated K.This option is very arbitrary, and in most cases it does not inhibit any recharge in porous and wet media (i.e.aquifers), as the recharge values are much smaller than the hydraulic conductivity.Third, we assume that preferential flow paths through soil can play a much larger role than the decreased hydraulic conductivity in most unconsolidated aquifers, as well as in most rock types (through faults and cracks).Therefore, and for the sake of simplicity of the NGRM model, it was chosen not to incorporate any calculations of unsaturated flow.
The rainfall recharge estimated by the NGRM model is a 'potential recharge', i.e., the recharge that would occur if the groundwater table is deep enough.Shallow groundwater tables will result in partial rejection of rainfall recharge, and a larger component might go to runoff or evaporation in these areas, which at this stage cannot be modelled by the NGRM model.Therefore, it is recommended to couple the NGRM with groundwater models, so that recharge can be corrected in areas such as wetlands or springs.This coupling process can then indicate areas where potential recharge can be corrected to actual recharge.
Finally, the NGRM (as well as the WaterGAP) model does not calculate the recharge of rivers to groundwater (except for rainfall on dry riverbeds).In some areas with large braided gravel river systems, losing rivers can recharge groundwater substantially, especially when streamflow is high and groundwater level is low (e.g., after heavy rainfall following a long dry period in autumn).
Recent research on nation-wide work on losing and gaining rivers is reported by Yang et al. (2015).
Therefore, future development of the NGRM is recommended to use results of that work.

475
For rainfall, the largest model input data component, NGRM application at the local scale reveals bias rather than uncertainty.For example, the national (VCS) rainfall data in the Waipa River catchment of the Waikato region, New Zealand, had to be increased by 15% by Rawlinson et al. (2015) in order to make it fit with the values of three independent models.Systematic errors or biases in model input components propagate differently from uncertainty as calculated by the NGRM model 480 equations.For example, a consistent bias in monthly rainfall, where there is under-catch, i.e., real rainfall is higher than measured rainfall, would lead to a different uncertainty of rainfall recharge model estimates.
Although ET is generally smaller than rainfall, the estimation of PET and AET also has considerable uncertainty.Westerhoff (2015) demonstrated that uncertainty in daily Penman PET in New

485
Zealand can be 10 -40%.For this study, it was assumed that the uncertainty of AET decreases with the AET/PET ratio.Because AET depends on more parameters than PET, e.g., soil moisture and vegetation health, the effect of the uncertainty of those parameters is unknown.This is one of the reasons that AET is often estimated from PET at the local scale, because soil parameters are considered to be better known at that scale.Furthermore, satellite-derived AET is known to contain errors  2006) (Figure 13), it does not necessarily mean it is better in every region.
More detailed research of canopy interception of the rainfall, and LAI, will improve model estimates.However, since the uncertainty of interception was assumed small in comparison with rainfall 495 and ET this will not be discussed here.
Comparison of the NGRM with the ECAN model in the Waimakariri catchment (see Results) showed that land surface recharge of the ECAN model is much higher in the river.This is not caused by a shortcoming of the model equation (as the ECAN model also does not embed a river recharge model), but by the differences in available soil data.The ECAN recharge model uses a local, more 500 spatially detailed, soil data set, which is not available at the national scale.These local data account better for recharge in the very permeable river areas.

Uncertainty of the NGRM in mountainous areas
The NGRM model indicates that rainfall recharge can be larger in mountainous areas, e.g., the flanks of the Southern Alps on the South Island (Figure 3).Although the existence of large exploitable 505 aquifers is unlikely in mountainous areas, rainfall recharge is likely to occur in these areas.This recharge is mostly caused by a high rainfall climate: although the geology seems relatively impermeable and much runoff is expected, substantial amounts of rainfall can still permeate into the ground.For example, Sims et al. (2015) surmise that a maximum of 20% of rainfall could infiltrate schist bedrock in the Southern Alps, New Zealand.These amounts could recharge to small, perched, 510 aquifers to enter as baseflow at relatively short distance.This is likely to occur in the higher regions of Canterbury or most of the West Coast on the South Island.Alternatively, groundwater could recharge to (unknown zones of) deep groundwater.For example, Calmels et al. (2011) find deep groundwater, likely to come from the mountainous region in Taiwan, with a similar high-rainfall climate and terrain to New Zealand.Doyle et al. (2015) show that 45% of aquifer recharge originates 515 from mountain blocks in their study area in British Columbia, Canada.
The uncertainty of rainfall recharge is, however, largest within mountainous regions.This is because, for low hydraulic conductivity, the large uncertainty of K plays a significant role in the recharge model uncertainty (e.g., Figure 5, right).Additionally, low permeability and high rainfall can cause overestimation of recharge in monthly estimates: the model assumes that this falls evenly 520 over the month, but the monthly rainfall is realistically the result of a few high-rainfall events that lasted not more than days or shorter.Daily estimates would improve recharge estimation in mountainous regions, but this would also impact the speed and simplicity of the model.Finally, recharge can be underestimated due to under-representation of fractures and faults in the model.Fractured zones, common in mountainous areas, are currently only partly embedded in the estimation of K (as 525 the method of Gleeson et al. (2011) is partly based on calibration of multiple hydrological models in hard-rock formations.
Comparisons with locally calibrated models in this research (Waimakariri and Mataura) show that the NGRM model estimates a low mean recharge (when uncertainty is not taken into account).As known, calibrated models do not take into account recharge that occurs outside the model boundaries.

530
Recharge in the foothills or in mountainous area is thus not taken into account in these models, since they are not within the model boundary.If recharge outside the model boundaries would occur in reality, but models do not take this into account, then those models would be calibrated wrongly.

Uncertainty of the NGRM in irrigated areas
In irrigated areas, the soil storage receives an additional irrigated amount of water.This is only partly incorporated in the NGRM.The effect of irrigation and interception is taken into account 540 by the AET, since the independent satellite-derived signal picks up vegetation health.Use of an independent satellite-derived signal is thus advantageous: it means that the AET is calculated as higher in these areas as the vegetation health has increased.Other parts of the model cannot always cope well with irrigation.For example, if irrigation is not fully efficient, (i.e., the water drains to groundwater instead of feeding the crop), the excess water will recharge and create an unknown  This discussion highlights the need for further research on application of adjusted large-scale models and satellite data on the national or regional scale.Future research for improvements in the NGRM model is summarised as: -Model improvements on rainfall-runoff, river recharge, snow and snowmelt, soil heterogene-560 ity, hydraulic conductivity; -Better uncertainty assessment of national input data, such as rainfall and AET; -Incorporation of larger (catchment-based) model boundaries in future local and regional recharge studies; -Correction of recharge through coupling with groundwater models in wetlands and springs;

565
-The effect of mountain recharge to groundwater modelling in New Zealand; -The added value of satellite-derived AET in irrigated areas.
The NGRM model is inspired by international, global-scale, researchers.If the NGRM model would be applied by regional councils, there is a need for further improvement of model equations and uncertainty estimates of both the model as well as its input components, e.g., rainfall and AET.Future research should therefore be performed in a catchment-by-catchment analysis with the best available data from regional councils, i.e., regional streamflow data and water budgets per subcatchment.In data-sparse regions, the regional data could then be completed with national flow data and statistics (e.g., Woods et al., 2006;Booker and Woods, 2014)        Hydrol.Earth Syst.Sci. Discuss., doi:10.5194/hess-2016-410, 2016 Manuscript under review for journal Hydrol.Earth Syst.Sci.Published: 30 August 2016 c Author(s) 2016.CC-BY 3.0 License. 115 140 used in New Zealand studies, since: they seem to fit expected values and patterns in large parts of New Zealand; and the data already take into account vegetation characteristics.The above described monthly MOD16 Penman PET and MOD16 AET, covering the period January 2000 to December 2013, were used in this study (mean annual AET is shown in Figure 1b).They are called PET and AET onwards. 145 220 Hydrol.Earth Syst.Sci.Discuss., doi:10.5194/hess-2016-410,2016   Manuscript under review for journal Hydrol.Earth Syst.Sci.Published: 30 August 2016 c Author(s) 2016.CC-BY 3.0 License.f soil is a correction factor for permeability of soil [0 to 1]; f geol is a correction factor for the geology [0 to 1].
rainfall recharge estimates 320 Application of the NGRM leads to nation-wide estimates of rainfall recharge at monthly intervals at 1 km grid resolution (mean annual values shown in Figure 3).The time series of recharge and 11 Hydrol.Earth Syst.Sci.Discuss., doi:10.5194/hess-2016-410,2016 Manuscript under review for journal Hydrol.Earth Syst.Sci.Published: 30 August 2016 c Author(s) 2016.CC-BY 3.0 License.some model components for three randomly chosen locations are shown in Figure 4. Mean annual uncertainty is highest at places where rainfall is high (Figure 5, left).However, relative uncertainty, i.e., as a percentage of the mean annual recharge (Figure 5, right), is largest in areas where hydraulic 325 conductivity is low.Recharge uncertainty is thus affected by rainfall for its high volumes, but most sensitive to geology for its high uncertainty and significant role in the model equation.The total New Zealand average rainfall recharge for the period 2000 -2014 estimated from the NGRM model is 2500 m 3 /s, with a model uncertainty σ of 17% (

335
NGRM recharge estimates were compared to lysimeters in the Canterbury Plains, New Zealand (Figure6).These four dryland lysimeter stations, located in the plains ('Airport', 'Hororata', 'Lincoln', and 'Winchmore'), have measured rainfall recharge from 1999 (or earlier) onwards.Hong and White (2014) compared these lysimeter observations with two locally applied models (called SOILMOD/DRAIN and SMB-SMC) for the period 2000-2004.Mean annual NGRM rainfall recharge 340 estimates for the same period are equal to rainfall recharge observations at three lysimeter stations (Table5 and Figure 7): recharge values fall within the model uncertainty σ at the Airport, Lincoln and Winchmore locations; at Hororata NGRM estimates much lower recharge, which is mainly caused by one anomalous rainfall event in January 2002 (Figure8), where torrential rains occurred on the 1st of January, and heavy rain wreaked havoc and caused surface flooding along the east coast 345 of Canterbury from 11-13 January(NIWA, 2002a, b).The locally calibrated SOILMOD/DRAIN model handles this anomalous wet event better at Hororata than the NGRM and SMB-SMC model.NGRM rainfall recharge estimates were compared to a locally calibrated model in the WaimakaririCanterbury Water Management Strategy Zone, New Zealand (Figure9).The regional council, Environment Canterbury (ECAN), developed an advanced land surface recharge model with a grid 350 resolution of 200 m for this areaAlkhaier (2016).This 'ECAN' model was built with the MIKE-SHE hydrological platform DHI (2016) and includes a sensitivity analysis with the PEST package(Doherty, 2016).The ECAN model gives mean recharge for the time period 1972 -2015 for three different model scenarios: a minimum, average, and maximum land surface recharge scenario (Figure 10).The 2000 -2013 mean annual recharge (196 ± 27 mm/yr) agrees well with the ECAN 355 model's minimum (195 mm/yr) scenario (

370
Recharge data were compared in the mid-Mataura catchment, located in the Southland Region, New Zealand (Figure11).The area for this case study consists of five groundwater management zones.The groundwater flow model developed for the mid-Mataura catchment(Burberry et al., 2013) used recharge values that were estimated for polygons with the Rushton method(Rushton et al., 2006).Mean annual recharge estimates between July 2000 and June 2007 of the NGRM and 375 Rushton model in the five groundwater zones have been used for this comparison.Total mean annual NGRM recharge is lower than the Rushton estimate (131 ± 27 mm/yr and 215 mm/yr, respectively, see Table 405 Hydrol.EarthSyst.Sci.Discuss., doi:10.5194/hess-2016-410,2016   Manuscript under review for journal Hydrol.Earth Syst.Sci.Published: 30 August 2016 c Author(s) 2016.CC-BY 3.0 License.

4704. 2
Limitation of model input dataWe considered that model input data of the NGRM model is much finer, and probably better, than that of the large-scale WaterGAP model.However, application of the NGRM model on the local scale, as should every model, should embed a careful consideration of the quality of model input data.

Future
local recharge models should therefore consider taking into account at least the scale of the whole catchment, i.e. including the foothills and mountains.It is recommended to perform more 535 research on mountain recharge, its relation to deep groundwater, and its relation to fracture zones, and implications for model calibration and model time steps.Hydrol.Earth Syst.Sci.Discuss., doi:10.5194/hess-2016-410,2016 Manuscript under review for journal Hydrol.Earth Syst.Sci.Published: 30 August 2016 c Author(s) 2016.CC-BY 3.0 License.
545bias in the monthly soil storage of the NGRM model.If water abstraction for irrigation comes from groundwater, the long-term effect of this excess irrigation recharge (a positive flux to recharge) will balance the water abstraction (a negative flux to recharge).However, if irrigation comes from surface water, this could impact both the monthly and the long term estimates of the NGRM recharge estimation.If national irrigation data become readily available, it is recommended to add these to the 550 model equation.However, estimating exact volumes of irrigation is not straightforward and could result in other, much larger, bias of the model estimates.Use of the satellite-derived AET could help in better estimates of irrigated and non-irrigated areas.Furthermore, use of the MOD16 PET in a ratio of AET to PET could give a better, measured, indication of soil water deficit, which could be fed back into the NGRM model.

Figure 10 .
Figure 10.Mean annual recharge in Waimakariri catchment management strategy zone.Left: NGRM recharge.Right: Recharge as evaluated by the ECAN land surface recharge model.

Figure 12 .
Figure 12.Mean annual recharge and its spatial distribution in the mid-Mataura catchment model area, according to the NGRM (left) and the Rushton (right) models.Cumulative probability in the bottom figures is shown in red.

Figure 13 .
Figure 13.Mean annual AET derived from lysimeters in the Canterbury Plains, New Zealand, compared to MOD16 AET and AET from Woods et al. (2006).

Table 4
).Total rainfall recharge values for the North Island is 1334 m 3 /s (a mean of 370 mm/yr, σ 15%); for the South Island it 330 is 1166 m 3 /s (a mean of 243 mm/yr, σ 18%).The national estimate falls within the 2835 m 3 /s (σ 10%) recharge estimate of

Table 6
rainfall recharge was in between approximately 0.25 and 0.45 of rainfall for the different scenarios.The spatial distributions of the mean annual recharge, quantified as a standard deviation, are 360 similar for the NGRM and the ECAN models (110 mm/yr for NGRM compared to approximately 125 mm/yr for the three ECAN scenarios).However, there are visual spatial differences between the NGRM and the ECAN recharge models, of which the clearest difference is in river areas, where ECAN recharge is high and NGRM recharge is low (see discussion section on river recharge).It seems more likely that the ECAN model is more realistic than the NGRM in these areas, since the ). Comparing recharge for the two different time periods only creates a maximum estimation error of approximately 12 mm/yr, since mean (VCS) rainfall for the model area for 1972 -2013 was 798 mm/yr; for 2000 -2013 it was 763 mm/yr; Hydrol.Earth Syst.Sci.Discuss., doi:10.5194/hess-2016-410,2016 Manuscript under review for journal Hydrol.Earth Syst.Sci.Published: 30 August 2016 c Author(s) 2016.CC-BY 3.0 License.and 365 ECAN model is calibrated.However, both models lead to the same mean annual recharge.To investigate which of the two models is best in which area requires more in-depth research that takes into account factors as: recharge in areas out of the model boundary (see discussion), groundwater table depths; spring locations, and baseflow separation methods at multiple gauging locations in the rivers.
Rushton recharge estimates exclude areas outside the model boundaries that are still in the catchment.This could have implications for recharge occurring in areas in the catchment, but outside the model boundary (see This is because of many reasons, of which some are explained here.First, the NGRM model has model limitations, mostly due to its simplifications.Although NGRM matches case studies well in this paper, the model may be too simple in other (sub-)catchments.These model limitations are described below in the topic 'Model equation limitations'.Second, current national input datasets might not Woods, 2014;Tait et al., 2006;Woods et al., 2006)hus include an uncertainty estimate of the Rushton recharge, and should include the whole (sub-)catchment and not only the model boundary.valuablenewadditiontoexisting national datasets of terrestrial water cycle variables (Booker and 390Woods, 2014;Tait et al., 2006;Woods et al., 2006).For example, national water budgets based on long-term means can be aided by: estimates of baseflow using the long-term mean rainfall recharge; estimates of AET that include irrigation and interception; estimates of unrouted quickflow using the Hydrol.Earth Syst.Sci.Discuss., doi:10.5194/hess-2016-410,2016Manuscript under review for journal Hydrol.Earth Syst.Sci.Published: 30 August 2016 c Author(s) 2016.CC-BY 3.0 License.surfacerunoff that is stored by the model; and estimates of rainfall that are corrected for terrain slope, interception, and snow.Theoretically, a water budget can be made using all NGRM data.However, 395 we acknowledge that implementation of a national water budget requires a much larger effort.400 be good enough for application at the local scale.This is discussed in the topic 'Limitations and advantages of nation-wide model input data'.Third, considerable uncertainty remains in data-sparse regions such as irrigated and mountainous areas, which is discussed in the topics 'Uncertainty of NGRM in mountainous areas' and 'Uncertainty of the NGRM in irrigated areas'.Recommendations are compiled and discussed in a final topic 'Future Research'.

Table 1 .
Gleeson et al. (2011)h organisations.This recommended research can be best applied in a collaborative environment, with regional 575 councils, national research organisations, and international researchers.In this context, all the mentioned topics of research lead to one overall recommendation, i.e., more and better collaboration between international research of large-scale models and data with national and regional stakeholders in the research fields of groundwater and surface water.International and national research funds should further stimulate this collaborative approach.Although the NGRM model is uncalibrated, its recharge estimates compare well to most local and regional lysimeter data and recharge models.From the case study comparisons it is concluded that 595 the nation-wide rainfall recharge model gives a valuable initial estimate when applied at the local or regional scale, and can thus also be used in areas as a valuable initial estimate in data-sparse areas.New Zealand hydrolithologies and their intrinsic permeability κ, including standard deviation σκ, afterGleeson et al. (2011)and improved for the New Zealand context.F-g.= fine-grained; p-s.= poorly sorted; c-g.= coarse-grained; uncons.= unconsolidated; sed.= sedimentary.
580 5 Conclusions This paper developed an approach to use large-scale, global, models and satellite data to estimate rainfall recharge at the national scale, i.e., across New Zealand, and smaller.The model, NGRM, uses an adjustment of the WaterGAP model, and MODIS derived ET and vegetation data, and furthermore the available nation-wide datasets on rainfall, elevation, soil and geology.The NGRM 585 estimates 1 km x 1 km monthly nation-wide rainfall recharge from January 2000 to December 2013.A valuable addition to the recharge estimation is the model uncertainty estimate (not typically output from previous rainfall recharge models in New Zealand), based on variance and covariance analyses, and model sensitivity of input components in the model environment.As opposed to the WaterGAP model, the NGRM model utilises better national input data of rainfall, soil, elevation, geology.The 590 estimated total New Zealand average recharge of the NGRM model results compiles to approximately 2,500 m 3 /s, or 298 mm/yr, with a model uncertainty of 17%, and is similar to the WaterGAP model.mostsensitive to rainfall in areas where recharge is high, but that uncertainty in hydraulic conductivity plays an important role in areas where recharge is impeded by geology.Future research topics for application of large-scale models at the smaller scale are recommended to focus on collabora-

Table 4 .
2000-2014 mean rainfall recharge values of the NGRM compiled for New Zealand, the North Island and the South Island, including model uncertainty.

Table 5 .
Mean annual rainfall recharge for July 2000 -June 2004 for four Canterbury lysimeters stations and the NGRM model.All values are in mm/yr.