Incorporation of Satellite Data and Uncertainty in a Nationwide Groundwater Recharge Model in New Zealand

A nationwide model of groundwater recharge for New Zealand (NGRM), as described in this paper, demonstrated the benefits of satellite data and global models to improve the spatial definition of recharge and the estimation of recharge uncertainty. NGRM was inspired by the global-scale WaterGAP model but with the key development of rainfall recharge calculation on scales relevant to nationaland catchment-scale studies (i.e., a 1 km × 1 km cell size and a monthly timestep in the period 2000–2014) provided by satellite data (i.e., MODIS-derived evapotranspiration, AET and vegetation) in combination with national datasets of rainfall, elevation, soil and geology. The resulting nationwide model calculates groundwater recharge estimates, including their uncertainty, consistent across the country, which makes the model unique compared to all other New Zealand estimates targeted towards groundwater recharge. At the national scale, NGRM estimated an average recharge of 2500 m3/s, or 298 mm/year, with a model uncertainty of 17%. Those results were similar to the WaterGAP model, but the improved input data resulted in better spatial characteristics of recharge estimates. Multiple uncertainty analyses led to these main conclusions: the NGRM model could give valuable initial estimates in data-sparse areas, since it compared well to most ground-observed lysimeter data and local recharge models; and the nationwide input data of rainfall and geology caused the largest uncertainty in the model equation, which revealed that the satellite data could improve spatial characteristics without significantly increasing the uncertainty. Clearly the increasing volume and availability of large-scale satellite data is creating more opportunities for the application of national-scale models at the catchment, and smaller, scales. This should result in improved utility of these models including provision of initial estimates in data-sparse areas. Topics for future collaborative research associated with the NGRM model include: improvement of rainfall-runoff models, establishment of snowmelt and river recharge modules, further improvement of estimates of rainfall and AET, and satellite-derived AET in irrigated areas. Importantly, the quantification of uncertainty, which should be associated with all future models, should give further impetus to field measurements of rainfall recharge for the purpose of model calibration.


Introduction
Satellite data have been shown to be important for global-scale assessments of water cycle variables, such as rainfall (GPM: [1]), evapotranspiration (ET) [2,3] and groundwater [4].Global-scale hydrological models are increasingly used to characterise Earth's groundwater [5][6][7], its fluxes [8,9] and its depletion [10].Applications of these data and models will become more important with predictions for more climate extremes (i.e., droughts [11] and floods [12]) and consequent pressures on agricultural production and hazard response.
Global-scale satellite data and hydrological models are mostly spatially coarse (typically with grid resolutions of hundreds of square kilometres) and therefore lack the reliability for application at smaller scales.High-resolution data are often needed to cover the diversity of landscapes and climates, such as occurring in countries like New Zealand, where mean annual rainfall can differ from approximately 0.5 to 10 m over relatively short distances and groundwater fluxes are highly variant over the nation.
Studies in spatially diverse climates have shown applicability of higher-resolution satellite data in national or catchment-scale parametrisation of ET.MODIS-based ET (with a resolution of approximately 1 km 2 ) was shown to be useful at the national scale of New Zealand [13] and Taiwan [14], and at the catchment scale, e.g., in Chili [15] and India [16].A recent study in Greece also showed that MODIS based ET showed comparable results to modelled AET [17].The well-known Landsat and Sentinel-2 satellite missions show applicability for more spatially refined (i.e., 100 m 2 grid resolution) ET estimates [18,19].
High-resolution satellite ET data might increase capability of global-scale models (e.g., groundwater recharge) towards application at the national or catchment scale.Because that data-model combination provides large coverage with high spatial resolution, potential benefits include: increasing the spatial resolution of global-scale models; filling the gaps in areas where ground observations are sparse; and providing consistency to better address trans-boundary issues.Trans-boundary issues can occur if water allocation policies are not compatible across boundaries because neighbouring regions use differing data standards and models to assess allocation.The incorporation of higher resolution satellite data into global-scale models should improve understanding of global-scale model behaviour, but might also better assess uncertainty of local-scale models: as most models are calibrated in their local environment, not much is known about their uncertainty, especially near or beyond the boundaries of study areas that are smaller than a whole catchment.
In New Zealand, regional councils are responsible for policy on water allocation.Regional boundaries generally follow catchment boundaries.Groundwater recharge is often a key input to set allocation limits.National guideline limits for groundwater allocation with respect to recharge exist [20,21], but policy and the subsequent recharge assessments are managed at the regional scale.Development of a high-quality national assessment of groundwater resources, including recharge, is challenging.This is because recharge models are inconsistent across regions, partly because of the geography of the country, which includes high mountain ranges and large areas of recent volcanic deposits.To date, New Zealand does not have this consistent nationwide characterisation of groundwater recharge.Therefore, calculation of national-scale groundwater recharge is a key challenge for water resource management.The lack of such a nationwide approach hinders further scientific advances in national hydrological modelling.For example, the New Zealand national surface water model [22] could benefit from better information of groundwater fluxes, such as recharge.
Global recharge models are equipped to cover large areas and could thus fulfill the need for national recharge information.Incorporation of higher-resolution data into these models, e.g., satellite data or national datasets, might further improve the spatial characterisation of those global models.Incorporation of global models with high-resolution data might better the effect of spatial coarseness of those global models.
Uncertainty has not typically been assessed by groundwater recharge models in New Zealand.Most recharge models in New Zealand are some form of soil water balance [23][24][25], where differences can be caused by their different model equations (e.g., derivations of actual ET and rainfall-runoff ratios) or their input data.The impact of recharge uncertainty was demonstrated by White et al. [25], who used three groundwater recharge models for the Central Plains of Canterbury, New Zealand.They showed that estimated groundwater use as a percentage of recharge was highly model-dependent, and because of that dependency varied in between 63% and 80% in a relatively dry year.Another study in New Zealand [26] showed that the difference in three groundwater recharge models applied in one region was most likely caused by the inconsistency of input data and not by the differences in model equations.Groundwater recharge models in New Zealand would therefore highly benefit from a better uncertainty assessment.
This paper describes the calculation of a nationwide model of groundwater recharge from rainfall for New Zealand that incorporates MODIS satellite data of ET and vegetation alongside other national input datasets and is inspired by the global-scale WaterGAP recharge model [8].The resulting nationwide model calculates groundwater recharge estimates, including their uncertainty, consistent across the country.The principal aim of this research is to assess the applicability of satellite data to downscale global-scale models, i.e., to national and catchment scales.The results are evaluated with local models and ground observations and an analyses of the inherent model equation uncertainty.This paper concludes with a discussion on future research topics, e.g., the advantage that satellite data brings to estimate irrigated areas; recommendations for improved model assumptions; and the need for collaborative research on the topic of groundwater recharge.

General Description
The model in this study is called the Nationwide Groundwater Recharge Model (NGRM).The NGRM has a 1 km × 1 km model grid that covers most of New Zealand, in the national New Zealand Transverse Mercator (NZTM) projection.The NGRM applies monthly MODIS-derived satellite data of evapotranspiration, i.e., MOD16 [3,27].Earlier studies calculated uncertainty of MOD16 data by taking into account the uncertainty of both ground observations and satellite data in the Penman and Penman-Monteith equations [13].That uncertainty of MOD16 data was used as one of the inputs for the model uncertainty calculation.Additionally, satellite data of vegetation (MODIS C5 LAI/FPAR [28,29], Figure 1) were used to calculate interception.Satellite data are just one of many input in the NGRM, which also applies data from national databases of precipitation, topography, soil parameters, and geology (Figure 2).A detailed description of input components is in Appendix A.  The NGRM calculates rainfall recharge to groundwater in monthly time steps, currently covering the period of January 2000 to December 2014.The model uses a simple soil water balance: 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 [8], USGS-SWB [24], Rushton [23].The 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 recharge correction factors, 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 NGRM is uncalibrated.However, case study comparisons described within this article have been used to improve the nationwide model equations such that the model is considered to provide suitable estimates in these case studies.

Rainfall Recharge Model Equation
In monthly time steps, rainfall recharge R (in mm) was calculated for time step i as: where: If R i < 0, it is set to zero and the resulting storage deficit S i < S i−1 is used for the calculation of R i+1 .Appendix B gives a detailed description of pre-processing, model assumptions, and the assignment of correction factors of slope, soil and geology.

Uncertainty of Rainfall Recharge
From here onwards, 'top-down' uncertainty is defined as the NGRM model equation uncertainty; 'bottom-up' uncertainty is defined as the comparison of NGRM modelled output to ground observations or other (local) models.
Top-down uncertainty was calculated by error propagation of Equation (1), i.e., propagating the variance and covariance of all input components [30]: 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.Appendix C describes the individual components of the variance-covariance matrix.
Top-down model uncertainties were compiled for further analyses.First, mean monthly and mean annual values for each model cell were compiled for comparison with ground-observed data and other local models.Second, NGRM recharge values were binned into simple histograms.Third, NGRM recharge values were binned into 100 bins in between the minimum and maximum recharge, where each bin was plotted against mean and standard deviation of top-down uncertainty per bin.Subsequently, individual model input components (P i , f slope , AET, S, f soil and f geol ) were binned and plotted against top-down uncertainty.

National Rainfall Recharge Estimates
Application of the NGRM leads to nationwide estimates of rainfall recharge at monthly intervals at 1 km × 1 km grid resolution (mean annual values shown in Figure 3).The time series of recharge and 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).Uncertainty of geology plays a minor role in large aquifers, where the geology is permeable, because it does not constrain recharge in those areas.However, relative uncertainty, i.e., as a percentage of the mean annual recharge (Figure 5, right), is largest in areas where hydraulic conductivity is low, which is where the model equation is most sensitive to the high uncertainty of hydraulic conductivity (see Section 4 and Appendix B).Recharge uncertainty is thus affected by rainfall for its high volumes, but also sensitive to geology for its high uncertainty and significant role in the model equation.Jan-00 Jan-01 Jan-02 Jan-03 Jan-04 Jan-05 Jan-06 Jan-07 Jan-08 Jan-09 Jan-10 Jan-11 Jan-12 Jan-13 Jan-14 Jan-00 Jan-01 Jan-02 Jan-03 Jan-04 Jan-05 Jan-06 Jan-07 Jan-08 Jan-09 Jan-10 Jan-11 Jan-12 Jan-13 Jan-14 Jan-00 Jan-01 Jan-02 Jan-03 Jan-04 Jan-05 Jan-06 Jan-07 Jan-08 Jan-09 Jan-10 Jan-11 Jan-12 Jan-13 Jan-14   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% (Table 1).Total rainfall recharge for the North Island is 1334 m 3 /s (a mean of 370 mm/year, σ 15%); for the South Island it is 1166 m 3 /s (a mean of 243 mm/year, σ 18%).The national estimate falls within the 2835 m 3 /s (σ 10%) recharge estimate of the WaterGAP model [8].Although that estimate is for a different period (1972)(1973)(1974)(1975)(1976)(1977)(1978)(1979)(1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990), rainfall in the two analysis periods is similar (1972 to 1990: 1881 mm/year; 2000-2014: 1839 mm/year).The resulting recharge data of this study are available for the scientific community (NETCDF-CF1.6format), through an open data licence (CC BY-NC 4.0: https://creativecommons.org/ licenses/by-nc/4.0/),by contacting the authors.

Analyses of Binned Recharge and Top-Down Uncertainty
The mean recharge uncertainty has most values in between 0 and 10 mm/day, although maximum values occur up to 35 mm/day (Figure 6, left).Recharge uncertainty, relative to recharge, is highest for low values of recharge (Figure 6, right).Mean uncertainty for very low recharge can also be higher than 100% of the recharge value: although not plotted in Figure 6 (right), 26% of all uncertainty values are higher than their recharge (i.e., higher than 100%); and the mean uncertainty in between 0 and 1 mm/day recharge is over 200% of recharge.With increasing recharge, uncertainty diminishes to approximately 13% of recharge.Figure 7 shows the behaviour of mean recharge and mean uncertainty per input component.In general, rainfall and geology play a dominant role on both recharge and uncertainty.This confirms the findings from Figure 5: in general, recharge according to the NGRM model is most sensitive to the underlying geology and rainfall (Figure 7, left axes).This analyses also shows that the NGRM model equation of recharge is most sensitive to rainfall when rainfall is less than 10 mm/day.Although other input components (AET, soil permeability, PAW and terrain slope) play a less dominant role, they still influence the recharge, e.g., for AET in between 0 and 2 mm/day mean recharge still shows considerable variation (Figure 7, top right).Recharge uncertainty σ R is also mostly depending on the geology and rainfall (Figure 7, right axes).The behaviour of σ R with increasing rainfall is different than that of recharge.Mean recharge uncertainty is higher than mean recharge over the whole range of f geol , i.e., when geology limits recharge, uncertainty becomes significantly higher, also confirming earlier conclusions from Figure 5 (right).

Bottom-Up Uncertainty: Comparison of NGRM Recharge with Published Case Studies
NGRM recharge estimates were compared to lysimeters in the Canterbury Plains, New Zealand (Figure 8).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 [31] compared these lysimeter observations with two locally applied models (called SOILMOD/DRAIN and SMB-SMC and running at a daily time step) for the period 2000-2004.Mean annual NGRM rainfall recharge estimates for the same period are equal to rainfall recharge observations at three lysimeter stations (Table 2 and Figure 9): recharge values fall within the top-down 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 (Figure 10), where torrential rains occurred on the 1st of January, and heavy rain wreaked havoc and caused surface flooding along the east coast of Canterbury from 11-13 January [32,33].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 Waimakariri Canterbury Water Management Strategy Zone, New Zealand (Figure 11).The regional council, Environment Canterbury (ECAN), developed an advanced land surface recharge model with a grid resolution of 200 m for this area [35].This 'ECAN' model was built with the MIKE-SHE hydrological platform [36].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 12).The 2000-2013 mean annual recharge (196 ± 27 mm/year) agrees well with the ECAN model's minimum (195 mm/year) scenario (Table 3).Comparing recharge for the two different time periods only creates a maximum estimation error of approximately 12 mm/year, since mean (VCS) rainfall for the model area for 1972-2013 was 798 mm/year; for 2000-2013 it was 763 mm/year; and 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 similar for the NGRM and the ECAN models (110 mm/year for NGRM compared to approximately 125 mm/year 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.It seems likely that the ECAN model is more realistic than the NGRM in these areas, since the ECAN model uses local input data on soil.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, groundwater table depths; spring locations, and baseflow separation methods at multiple gauging locations in the rivers.Recharge data were compared in the mid-Mataura catchment, located in the Southland Region, New Zealand (Figure 13).The area for this case study consists of five groundwater management zones.The groundwater flow model developed for the mid-Mataura catchment [37] used recharge values that were estimated for polygons with the Rushton method [23].Mean annual recharge estimates between July 2000 and June 2007 of the NGRM and 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/year and 215 mm/year, respectively, see Table 4).This is partly caused by the difference in rainfall used for both models (809 mm/year for the NGRM model vs. 903 mm/year for the Rushton model).However, the ratio of recharge and rainfall of NGRM (0.16 ± 0.03) is still lower than the Rushton model (0.24).If Rushton and NGRM uncertainties would be similar, then recharge estimates NGRM would be 2σ lower than the Rushton output.However, the Rushton polygons have an unknown but potentially large uncertainty, since there are so few polygons compared to the NGRM model cells (Figure 14).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 discussion).Further local comparison should thus include an uncertainty estimate of the Rushton recharge, and should include the whole (sub-)catchment and not only the model boundary.

Satellite Data to Better Constrain 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 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 bias 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 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.

Limitations of the Model Equation and the Model Input Data
NGRM rainfall recharge estimates, including their underlying input satellite data of AET and vegetation, are a valuable new addition to existing national datasets of terrestrial water cycle variables [38][39][40].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 satellite-derived AET that include irrigation and interception; estimates of unrouted quickflow using the surface runoff 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, we acknowledge that implementation of a national water budget requires a much larger effort.This is because of many reasons, of which most relate to either model equation limitations or model input data.
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.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 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 on, e.g., slope-runoff relations, snow correction, hydraulic conductivity and river recharge (see Appendix A).Although case study comparisons have been used to improve the nationwide model equations such that the model is considered to provide suitable estimates in these case studies, uncertainty in explored areas is still relatively unknown.For example, considerable uncertainty remains in mountainous areas.Despite mountain geology being relatively impermeable, substantial amounts of rainfall can still permeate into the ground in high-rainfall climates (e.g., Southern Alps, New Zealand [41], Taiwan [42], British Columbia, Canada [43]).However, the uncertainty of rainfall recharge in mountains is large, because of the large uncertainty of hydraulic conductivity estimates in the NGRM model equation (see Results) and the often unknown location and permeability of fractures and faults.Comparisons with locally calibrated recharge models in this research (Waimakariri and Mataura) show that the NGRM model estimates a low mean recharge (when uncertainty is not taken into account).These local models do not take into account recharge that occurs outside the model boundaries.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 could be applied or calibrated incorrectly.Local recharge models should therefore consider taking into account at least the scale of the whole catchment, including the foothills and mountains.
Current national input datasets might not be good enough for model application at the local scale.We considered that model input data of the NGRM model is much finer, and probably better, than that of the large-scale WaterGAP model.However, local application of the NGRM model should embed a careful consideration of the quality of model input data (as should every model).This paragraph discusses three important model input data components: rainfall, ET and soil.Rainfall: NGRM application at the local scale reveals bias in rainfall rather than uncertainty.For example, the national (VCS) rainfall data in a study in the Waipa River catchment of the Waikato region, New Zealand, had to be increased by 15% in order to make it fit with the values of three independent models [26].Systematic errors or biases in model input components propagate differently from uncertainty as calculated by the NGRM model equations.ET: although MOD16 AET (2000-2013) compares better to lysimeter-derived AET (1999-2011) in the Canterbury Plains of New Zealand than the commonly used average 1960-2006 AET [40] (Figure 15), an earlier study [13] demonstrated that model 'top-down' uncertainty in daily Penman PET in New Zealand can be 10-40% and estimation to AET brings more uncertainty because soil moisture and vegetation health are not always well known.Soil: 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 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.

Future Research
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:

•
The added value of satellite-derived AET in irrigated areas; • Model improvements on rainfall-runoff, river recharge, snow and snowmelt, soil heterogeneity, 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;

•
The effect of mountain recharge to groundwater modelling in New Zealand; The NGRM model is inspired by international, global-scale, research.However, application of the NGRM model at th regional scale requires further improvement of model equations, input data and uncertainty estimates, e.g., locally measured rainfall and ET.Future research should therefore include all available data, starting with the best available, i.e., data from regional councils.In data-sparse regions, the regional data could then be completed with nationwide data, including data from this study and national flow data and statistics (e.g., [38,40]).This recommended research can be best applied in a collaborative environment, with regional 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.

Conclusions
This paper developed an approach to use large-scale satellite data and global hydrological models to estimate rainfall recharge at the national scale, i.e., across all regions of New Zealand.The model, NGRM, is inspired by the global WaterGAP model, but uses higher resolution input data , i.e., MODIS derived ET and vegetation data and the available nationwide datasets on rainfall, elevation, soil and geology.The NGRM estimates 1 km × 1 km monthly nationwide rainfall recharge from January 2000 to December 2014.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.The estimated total New Zealand average recharge of the NGRM model results compiles to approximately 2500 m 3 /s, or 298 mm/year, with a model uncertainty of 17%.Those results are similar to the WaterGAP model, but the improved input data results in much better spatial characteristics of recharge estimates.Total rainfall recharge for the North Island is 1334 m 3 /s (a mean of 370 mm/year, σ 15%); for the South Island it is 1166 m 3 /s (a mean of 243 mm/year, σ 18%).
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 the nationwide 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.Local applications might require the model to be calibrated and, as with any model, it is therefore recommended to carefully consider the NGRM model limitations for local application.
This research also provides improved insights into the uncertainty of rainfall recharge models, including the role of recharge model input components.It shows that recharge is most sensitive 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.Thus, the satellite data of ET and vegetation improve the spatial characteristics of the recharge model without adding significant uncertainty.This research also shows that the combined top-down and bottom-up approach gives a more insightful description of uncertainty and underlines the importance of the availability of ground observations.Future research topics for application of large-scale models at the smaller scale are recommended to focus on collaborative research between international and national researchers and regional water managers, thus covering all spatial scales.Topics for future research include, amongst others: the added value of satellite-derived AET in irrigated areas; improvement of estimates of rainfall and evapotranspiration; and the impact of fracture zones on rainfall recharge within mountainous areas.publish in open access have been covered by the SAC programme.This project has received co-funding from the European Union's Seventh Programme for research technological development and demonstration under grant agreement No. 603608, eartH2Observe.We furthermore would like to thank: Andrew Tait from NIWA for generously sharing gridded VCS data; Taejin Park, PhD student at Boston University, Department of Earth and Environment, Cliveg Research Group, for sharing the LAI dataset; Fouad Alkhaier (Environment Canterbury) for sharing recharge data in the Waimakariri case study; and Catherine Moore (GNS Science) for providing data of the mid-Mataura catchment and inspiration for the uncertainty analyses.
profile [55].PAW data are available from the LRIS portal as an ArcGIS shapefile, with each soil class having a minimum and maximum value (Table A3).Modal values were calculated according to [54] and shown in Figure 2f.For areas where PAW was unknown, the mean of all PAW modal values was assumed.

Appendix B. Rainfall Recharge Model Equation (Detailed)
In monthly time steps, rainfall recharge R (in mm) was calculated for time step i as: where: If R i < 0 (Equation (1)), then R i is set to zero and the resulting storage deficit S i < S i−1 is used for the calculation of R 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 due to slope: where the error function (erf) fits the empirical slope-runoff relation used by the WaterGAP model [8].
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., [23,24]).Also, the NGRM model prefers the f slope to directly affect rainfall, and not affect AET and S.
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) [25,58] 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 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 A2).For f soil , it was assumed that soils lower than a value of 0.15 (equivalent to a soil permeability of less than 4 mm/h, [54]) reject 75% of the rainfall recharge to surface runoff; all other soils accept all recharge.The model assumes an isotropic permeability, i.e., 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 rainfall recharge.The factor f geol was then calculated as the ratio (in between 0 and 1) of potential rainfall recharge and hydraulic conductivity.

Appendix B.1. Model Equation Limitations
The slope-runoff relation is based on a sparse dataset of empirical values of the WaterGAP recharge model [8], 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 number (CN) approach relating soil type and land cover to runoff [59].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 [23,24] without taking into account terrain slope.The NGRM prefers to use the simple relation to terrain slope, because of 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' [60] 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 [61,62], leaving less erodible 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 [63].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 2014 now takes less than an 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 × 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 [60,64].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.Recent research also acknowledged this model limitation [65].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 [66].Therefore, future development of the NGRM is recommended to use results of that work.all input components (P, 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 rain gauges) to 15% in steep-sloped regions (without rain gauges).This is in line with earlier findings of maximum uncertainty of 15% [39]; 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., [67]) and averaging of the terrain model to the model grid;

•
The uncertainty in daily Penman PET in New Zealand can be 10-40% [13] 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 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.Uncertainty in K can be very high: even higher than the actual value of K [5,68].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).
Rescaling the covariance matrix to the values of the uncertainty in the input components was done as follows:

Figure 2 .
Figure 2. Overview of NGRM model input components of (a) precipitation (compiled to mean annual); (b) AET (compiled to mean annual values); (c) elevation; (d) hydraulic conductivity K; (e) soil permeability (here shown as a ratio, see Appendix A); and (f) PAW.

Figure 3 .
Figure 3. Nation-wide rainfall recharge from the NGRM compiled to mm/year.

Figure 5 .
Figure 5. Uncertainty σ rech of mean annual rainfall recharge estimated with the NGRM model as mm/year (left) and as a percentage of total recharge (right).

Figure 6 .
Figure 6.Left: Histogram distribution of all NGRM recharge uncertainty values, converted to mm/day; Right: the behaviour of mean recharge uncertainty over recharge, plotted over two axes: in mm/day (left axis, blue), including the standard deviation 1σ (grey band); and as a percentage of recharge (right axis, brown, 1σ not plotted for better visual convenience).

Figure 7 .
Figure 7. NGRM model input components plotted against recharge (left axis) and top-down uncertainty (right axis).

Figure 11 .
Figure 11.Model area in the Waimakariri catchment management strategy zone, Canterbury, New Zealand.

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

Figure 14 .
Figure 14.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 15 .
Figure 15.Mean annual AET derived from lysimeters in the Canterbury Plains, New Zealand, compared to MOD16 AET and the commonly used long-term average (1960-2006) AET from [40].

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

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

Table 3 .
Mean annual recharge for the Waimakariri catchment management zone model area (mm/year).
xy is the variance from the covariance analysis; • σ x,new and σ y,new are the errors in the input components; • σ x and σ y are the standard deviations from the covariance analysis;