Calibration of a transient transport model to tritium data in streams and simulation of groundwater ages in the western Lake Taupo catchment, New Zealand

Here we present a general approach of calibrating transient transport models to tritium concentrations in river waters developed for the MT3DMS/MODFLOW model of the western Lake Taupo catchment, New Zealand. Tritium has a known pulse-shaped input to groundwater systems due to the bomb tritium in the early 1960s and, with its radioactive half-life of 12.32 yr, allows for the determination of the groundwater age. In the transport model, the tritium input (measured in rainfall) passes through the groundwater system, and the simulated tritium concentrations are matched to the measured tritium concentrations in the river and stream outlets for the Waihaha, Whanganui, Whareroa, Kuratau and Omori catchments from 2000–2007. For the Kuratau River, tritium was also measured between 1960 and 1970, which allowed us to fine-tune the transport model for the simulated bomb-peak tritium concentrations. In order to incorporate small surface water features in detail, an 80 m uniform grid cell size was selected in the steady-state MODFLOW model for the model area of 1072 km 2. The groundwater flow model was first calibrated to groundwater levels and stream baseflow observations. Then, the transient tritium transport MT3DMS model was matched to the measured tritium concentrations in streams and rivers, which are the natural discharge of the groundwater system. The tritium concentrations in the rivers and streams correspond to the residence time of the water in the groundwater system (groundwater age) and mixing of water with different age. The transport model output showed a good agreement with the measured tritium values. Finally, the tritium-calibrated MT3DMS model is applied to simulate groundwater ages, which are used to obtain groundwater age distributions with mean residence times (MRTs) in streams and rivers for the five catchments. The effect of regional and local hydrogeology on the simulated groundwater ages is investigated by demonstrating groundwater ages at five model cross-sections to better understand MRTs simulated with tritium-calibrated MT3DMS and lumped parameter models.


Introduction
There is still limited understanding of the dynamics of the groundwater component that transmits much of the water from rainfall to streams (Loague et al., 2005;Ebel et al., 2008). Better understanding of groundwater flow paths and travel times is essential to improve understanding of movement of contaminants (e.g. nitrate) through the groundwater system into surface waters such as streams and lakes (McMahon et al., 2010;Stichler et al., 2008). Tritium, a direct tracer of groundwater dynamics (groundwater age and mixture of groundwater with different ages), can be used to calibrate quantitative groundwater flow and tracer transport models and thus to greatly improve our conceptual understanding of baseflow generation. Simulating groundwater flow and tracer transport with distributed parameter models (DPMs) allows high level of complexity such as spatially distributed groundwater recharge, surface water networks, aquifer hydrogeology, and chemical input into groundwater.  Zuber et al. (2011) identified two approaches of calibration to measured groundwater tracer data: (1) use of groundwater ages or mean residence times (MRTs) estimated from measured concentrations with lumped parameter models (LPMs), and (2) use of measured tracer concentrations directly for calibration of DPMs. In this second approach, steady-state or transient groundwater flow models are calibrated to measured groundwater heads and baseflows by adjusting hydraulic conductivity, groundwater recharge and stream bed conductance in the model. Then, transport model results are matched to a groundwater age tracer concentrations to calibrate aquifer porosity and dispersivity in the model, and to make some additional adjustments to the previously calibrated groundwater flow model parameters (Zuber et al., 2011). In both approaches, a groundwater age tracer is required for the transport model calibration meeting the following conditions: (1) defined input concentrations to the groundwater system, (2) non-reaction with aquifer materials and dissolved chemicals in groundwater, and (3) accurate measurements at sampling points of monitoring bores and streams.

Published by Copernicus
Many groundwater tracers have been sampled in bores, and their applicability to characterize groundwater flow and transport in aquifers is well established. However, the choice of a calibration tracer is not obvious for streams. The gaseous tracers such as SF 6 , chlorofluorocarbons (CFCs), and 3 H-3 He, which have well-defined atmospheric concentrations, cannot be easily applied as calibration tracers for DPMs due to open atmosphere exchange in stream waters. The dissolved chemicals such as chloride, sulphate, and nitrate measured in streams lack input information at the groundwater table or have no time-dependent input, and can undergo chemical reactions in the aquifer and stream waters (McMahon et al., 2010).
Tritium seems to be the ideal conservative tracer for the calibration of DPMs in streams. Tritium input concentrations are well established globally from the Global Network of Isotopes in Precipitation, which is run by the International Atomic Energy Agency, including many stations in NZ, the South Pacific, and Australia. Tritium is part of the water molecule, and its isotopic concentration in groundwater is not altered by chemical reactions in the aquifer. This makes tritium a conservative tracer in groundwater systems. Because of its radioactive decay, tritium concentrations in groundwater are dependent on the travel times, even at current constant tritium input concentrations in the Southern Hemisphere and faded bomb tritium (Morgenstern et al., 2010;Stewart et al., 2012). A similar situation holds in the Northern Hemisphere, where tritium concentrations in groundwater are rapidly declining but have not yet returned to the natural tritium levels of the pre-bomb era.
To our knowledge, no previous work has been published on calibration of DPMs using tritium data in streams. Several studies have simulated tritium concentrations in monitoring wells using transport models, which incorporate advection, dispersion, and chemical reactions of contaminants in groundwater systems. Zuber et al. (2005) simulated tritium and SF 6 concentrations in five observation wells in the 200-m-deep Bogucice multi-aquifer system with a MODFLOW/MT3DMS finite-difference model in southern Poland. The study area of 172 km 2 was represented by 250 m by 250 m grid cells with five layers in the finite-difference grid using Visual MODFLOW (VMOD) graphical user interface (Zuber et al., 2005). Zoellmann et al. (2001) simulated tritium and SF 6 with MODFLOW/MT3DMS in nine monitoring wells installed in a basalt aquifer in Germany. In addition, transport models were used to simulate 18 O and 2 H concentrations - Krabebenhoft et al. (1990) identified groundwater and lake connectivity by simulating 18 O transport in the four-layer finite-difference model. Stichler et al. (2008) estimated Lake Leis contribution to pumping wells by simulating 18 O and 2 H concentrations in the transient FEFLOW model. Tritium concentrations in bores have been simulated with particle tracking MODPATH (Szabo et al., 1996;Pint et al., 2003;Boronina et al., 2005;Hunt et al., 2006;Walker et al., 2007;Fienen et al., 2009;Starn et al., 2010;Eberts et al., 2012), but no work with MODPATH has been carried out on streams. MODPATH is a post-processing program that uses a completely different mechanics of modelling tritium concentrations from MT3DMS (Pollock, 1994). In MODPATH, an advective component of groundwater flow is considered that does not simulate concentrations; tritium concentrations are obtained by taking a convolution integral of transit time distributions. MT3DMS, on the other hand, solves the transport equation in single-or dual-domain to obtain simulated concentrations.
In this study a general approach of DPM calibration using tritium data in streams is presented and the benefits of this approach are illustrated. We use the second calibration approach by Zuber et al. (2011) with measured tritium concentrations in river and stream waters for the finite-difference groundwater flow (MODFLOW) and transport (MT3DMS) model of the western Lake Taupo catchment (WLTC). The steady-state MODFLOW model was calibrated to groundwater levels and baseflow observations. The transient MT3DMS model was calibrated with the measured tritium values in streams and rivers for the Waihaha, Whanganui, Whareroa, Kuratau and Omori catchments. Then, the transport model calibrated to tritium data in streams was used to simulate groundwater ages for the WLTC enabling us to gain understanding of the hydrogeology and groundwater-stream interaction in the WLTC.

MODFLOW model setup
In the conceptual model for WLTC, the aquifer system receives recharge from rainfall and drains into Lake Taupo via streams or directly through the lake bed. In Fig. 1, Oruanui Ignimbrite, Whakamaru Ignimbrite, Pakaumanu Ignimbrite and Greywacke basement units (listed from the most to least conductive) were extracted from the existing 3-D geologic model for the WLTC area (White et al., 2009). The Whakamaru unit was subdivided into "a" and "b" areas that are separated by a hidden fault based on the geological map for the WLTC area. These units were represented by 16 grid layers with uniform thickness of 20 m, and each grid cell was assigned the hydrogeological properties of one unit. The model area of 1072 km 2 was represented by a finite-difference grid of 500 rows and 335 columns with a uniform grid cell size of 80 m to allow small surface water features to be included in the model albeit with an exaggerated width (Haitjema et al., 2010).
The top elevation of the model grid was obtained from the modified Lake Taupo bathymetry and digital terrain model (DTM) data, and processed to reduce higher elevations in adjacent cells that exceed a specified vertical difference of 10 m and to have slopes of no more than 12.5 %. This DTM adjustment allowed us to reduce the slope of the terrain in steep slope areas such as ridges and cliffs while maintaining a constant layer thickness in the entire model. In layer 1, the Lake Taupo water level of 357 m was incorporated as a constant head boundary, and surface water features such as rivers were incorporated as drain cells with uniform conductance of 40 m 2 day −1 (Fig. 1). The locations of the drain cells were obtained from the 1 : 50 000 topographic dataset based on vector data of the streams, and the assigned drain bottom elevation was 1 m below the top elevation of layer 1. Groundwater recharge was assigned to the topmost saturated model layer (McDonald and Harbaugh, 1988;Harbaugh et al., 2000). The total groundwater recharge values were obtained from the difference of average annual rainfall  and actual evapotranspiration data (Tait and Woods, 2006) between 1960 and 2006 and reclassified into 10 zones with a discrete recharge value from 400 to 1300 mm yr −1 with 100 mm yr −1 increments.
The model domain was divided into 49 unique nonoverlapping zones to include 30 river gauging stations using ZONEBUDGET (ZBUD), which is a post-processing water flow and mass budget utility for MODFLOW (Harbaugh, 1990). The simulated river and stream baseflows were extracted from the ZBUD output and routed from headwaters towards the lake to obtain cumulative baseflows at gauging stations. The steady-state groundwater flow model was calibrated to median values of stream baseflows at 30 gauging stations and of groundwater levels at 18 bores, which are screened in layers 1-2 except for the deepest borehole of 58 m screened in layer 3. As a result of the trial-and-error calibration to groundwater levels (R 2 = 0.992; normalized RMS = 3.7 %) and stream baseflows (R 2 = 0.9695), all input recharge values shown in Fig. 1 were multiplied by 0.88. The hydraulic conductivity values in the calibrated groundwater flow model are 1 m day −1 for the Oruanui unit, 0.2 m day −1 for the Whakamaru unit (W5a), 0.08 m day −1 for the Whakamaru (W5b) and Pakaumanu (P7) units and 0.01 m day −1 for the greywacke. These calibrated values are comparable to the values measured in the field. The results of the groundwater flow model are shown in Fig. 2. Simulated groundwater elevations range from 357 m above mean sea level at the Lake Taupo lakeshore to over 1000 m in the northern part of the model domain and are influenced by the surface water network. Groundwater divides occur at elevated terrain areas below the model top layer. The depth to water table in plan and cross-sectional views indicates potential thickness of the unsaturated zone.

Tritium input and measurements
Monthly concentrations of tritium in rainfall were aggregated to annual averages and assigned as annual recharge concentrations in MT3DMS from 1950 to 2012. Tritium concentrations are reported as the tritium ratio (TR) or tritium units (TU), and the tritium analysis details are given in Morgenstern and Taylor (2009). From 1960, the rainwater has been sampled monthly for tritium concentrations at the Kaitoke station near Wellington. The tritium concentrations between 1955 and 1960 are deduced from waters in the Hutt River (Morgenstern and Taylor, 2009) and prior to 1955 ambient tritium concentration of 1.9 TR. Tritium was also measured at several rainfall stations across New Zealand to establish the scaling factor for the various regions to the Kaitoke record according to latitudinal and altitudinal variations over New Zealand. The scale factor is 1.0 for the Taupo region (Morgenstern, 2007), and for the Lake Taupo climatic conditions the annual tritium concentrations are not significantly changed during the recharge process due to seasonal variation of infiltration (Morgenstern et al., 2010).
Tritium concentrations measured in stream and river waters of the Waihaha, Whanganui, Whareroa, Kuratau and Omori catchments were used as calibration targets for the MT3DMS model. The tritium sampling in streams was carried out during baseflow conditions, when stream baseflow was commonly considered to be sourced predominantly from groundwater contributions. Tritium data from 2001 and 2002 are obtained from Vant and Smith (2002) and from 2004 to 2007 from Morgenstern (2007). In addition, tritium data from 1964, 1966 and 1970, which is the time of the tritium peak in rainwater due to nuclear weapons testing in the early 1960s, are shown in Morgenstern and Taylor (2009) for the Kuratau River. These data allow us to resolve the bomb-peak ambigu-ity in age interpretation (Stewart et al., 2012), and to fine-tune the transport model.

MT3DMS model setup
In the MT3DMS model, tritium was specified as a single species with first-order decay of 12.32 yr (Morgenstern and Taylor, 2009). Both single-and dual-domain transport model settings were used for tritium transport simulations in fractured ignimbrite rocks using MT3DMS version 5.2 (Neville, 2006;Zheng and Wang, 1999). In the single-domain model, tritium concentrations occur in flowing groundwater. For the dual-domain model, stagnant groundwater is present with tritium concentrations that are exchanged with the mobile phase using the mass transfer coefficient in the MT3DMS model. The high and low values of the mass transfer rate make the dual-domain model equivalent to the single-domain model with combined porosity and the single-domain model with the porosity of the mobile domain, respectively (Zheng and Wang, 1999). The immobile domain was assigned uniform porosity of 0.001, and the first-order mass transfer coefficient between mobile and immobile domains was estimated to be 1.54 × 10 −5 based on the values of the sphere shape of the fractures and the molecular diffusivity and porosity of the immobile domain (Neville, 2006;Bear and Cheng, 2010). The effective porosity was used to represent the porosity of the mobile domain in the model and was obtained from the tritium calibration to be 0.1 [-] for the Oruanui unit (O3a), 0.035 for the Whakamaru unit a (W5a), 0.006 for the Whakamaru unit b (W5b) and the Pakaumanu unit (P7), and 0.001 for the greywacke unit (B2).
The longitudinal dispersivity was selected to be 10 m, and the ratios of horizontal to longitudinal dispersivity and of vertical to longitudinal dispersivity are 0.1 and 0.01, respectively (Zheng and Wang, 1999). An implicit GCG solver with upstream finite-difference advection was selected, and the time steps for the MT3DMS simulation were specified with a Courant number of 0.75, a step multiplier of 1.0, and a step size of 365 days (Zheng and Wang, 1999). The total time span simulated with MT3DMS was 62 yr to simulate tritium concentrations from 1950 to 2012. Annual tritium concentrations were assigned as recharge concentrations in the MT3DMS model using a MTH recharge concentration boundary condition file for VMOD that was prepared with a custom Python code. The background tritium concentration in precipitation was assigned as recharge concentration in the MT3DMS, and the tritium transport model was run for 200 yr to archive background tritium concentrations in the aquifer, which were used as initial tritium concentrations for the 62-yr MT3DMS run.
The tritium concentrations in the stream and river waters were obtained by extracting simulated tritium concentration values at drain cells in each of five river catchments using Python script. The Python script allowed us to interrogate the binary array output fluid (BGT) and concentration (UCN) files. For each drain cell, the amount of tritium was obtained by multiplying the simulated tritium concentration at that cell by the simulated groundwater outflow. The average tritium concentration at the stream and river catchment outflow was obtained by weighting the amount of simulated tritium for a catchment by the total simulated baseflow at that catchment.
In addition, tritium concentrations in river waters were obtained by using the fluid budget ZBUD and mass budget MT3MS Transport Observation Package (TOB). The MOD-FLOW drain cells in each catchment were selected for the MT3DMS TOB mass flux object (Zheng, 2006). The groundwater flow MODFLOW and transport MT3DMS models were run using "translate and run" in VMOD. The translated MT3DMS model files by VMOD are used to run the MT3DMS model with mass flux output file from command prompt. This was required due to inability of VMOD to include the TOB package for MT3DMS. The mass of tritium reported by the TOB mass flux object reports is converted to tritium concentration with the use of simulated baseflow for each of the five river catchments. However, the ZBUD and TOB approach is more cumbersome to set up for large catchments.

Simulated tritium concentrations
The measured and simulated tritium concentrations for the Kuratau River waters are shown in Fig. 3. The measured tritium concentrations in Kuratau River water are shown by black dots with vertical lines that indicate the one-σ measurement errors and are within 0.04-0.06 range after year 2000 (Morgenstern and Taylor, 2009). The simulated tritium output using the calibrated dual-domain MT3DMS model is shown by the solid black line (RMS = 1.62). The dualdomain model simulations matched the measured tritium time series data best in both 1960-1970 and 2000-2008. The single-domain results matched 1960-1970 measurements, but were too low for 2000-2009 measurements. The measured tritium values between 1960 and 1970 were essential to fine-tune the transport model to produce simulated tritium concentrations during the bomb peak. The importance of bomb-peak tritium data is demonstrated by the dashed and dashed-dotted lines. The dashed line shows the simulated tritium outputs for porosity values of 0.07 for the Whakamaru unit a (W5a), 0.012 for the Whakamaru unit b (W5b) and the Pakaumanu unit (P7) and is everywhere below the tritium outputs simulated with the calibrated model. The dasheddotted line shows the simulated tritium outputs for porosity values of 0.017 for the Whakamaru unit a (W5a), 0.003 for the Whakamaru unit b (W5b) and the Pakaumanu unit (P7), and is everywhere above the tritium outputs simulated with the calibrated model. These largest differences between the dashed and the dashed-dotted lines are seen at the bomb peak in Fig. 3.
The tritium outputs simulated with the MT3DMS transport model are shown in Fig. 4 for the Waihaha, Whanganui, Whareroa, Kuratau and Omori catchments. For all five catchments, the simulated tritium concentrations show similar patterns in 2001-2009 and the largest difference of simulated tritium concentrations at the bomb peak (Fig. 4). Simulated tritium outputs are highest for Omori River with simulated flow of 0.5 m 3 s −1 from a drainage area of 27.4 km 2 , and lowest for the Whareroa Stream with simulated flow of 1.4 m 3 s −1 from an area of 59.3 km 2 . These simulated tritium outputs are in agreement with the measured tritium concentrations for the Whareroa Stream with the lowest tritium concentrations, shown by blue diamonds, and the Omori River with the highest concentrations, shown by pink downward triangles. The Whanganui Stream and Kuratau River have similar simulated tritium outputs although their catchments have different measured tritium data and hydrogeologic settings. For example, the Whanganui Stream has simulated flow of 1.6 m 3 s −1 from an area of 65.3 km 2 and measured tritium concentrations shown by red squares. The Kuratau River has simulated flow of 5.6 m 3 s −1 from an area of 199.6 km 2 . The Waihaha River with simulated flow of 3.9 m 3 s −1 from an area of 155.4 km 2 has simulated tritium concentrations that fall between those of Whareroa Stream and Kuratau River and measured tritium concentrations shown by green upward triangles.
The influence of groundwater and transport model parameters on simulated tritium concentrations in rivers was investigated following Zuber et al. (2011). In all comparisons conducted, model parameters such as dispersion, aquifer porosity, hydraulic conductivity, conductance and rainfall recharge had the most influence on the simulated tritium concentration during the bomb peak. In advective dominated aquifer systems such as WLTC, the dispersion and immobile domain porosity have very little influence on the simulated tritium concentrations in rivers. For example, an order of magnitude increase of dispersion relative to the tritium-calibrated model resulted in an up to 4 % reduction of simulated tritium concentrations in Whareroa between 1960 and 1970. An order of magnitude decrease of immobile domain, which decreased the volume of stagnant water in fractures, resulted in an up to 4 % increase of the simulated tritium concentrations in Omori. A 25 % increase of mobile domain porosity values decreased simulated tritium concentrations by up to 8 % in Waihaha. This is due to slower groundwater flow causing increases in transit times compared to the calibrated model, and therefore tritium has more time to decay, especially at deep layers. Consequently, reduction of model layers produces an increase in simulated tritium concentrations due to shorter groundwater transit times in the aquifer sys-tem. In the groundwater flow model, a 25 % reduction was applied to groundwater recharge, conductance and hydraulic conductivity values. This reduction in groundwater recharge results in an up to 20 % decrease of the simulated tritium concentrations in Whareroa during the bomb peak. This is due to the smaller amount of tritium entering the aquifer system with groundwater recharge and lesser groundwater velocities that increase transit time. The reduction of the conductance and hydraulic conductivity values resulted in up to 3 % and up to 7 % decrease in Whareroa during the bomb peak, respectively. To understand these complex groundwater flow patterns, we conducted groundwater age simulation with the tritium-calibrated MT3DMS model.

Simulated groundwater age with tritium-calibrated model
Groundwater residence times at sampling points such as groundwater wells, springs and rivers have been determined from isotope measurements using lumped-parameter (box) models (LPMs) such as exponential, piston and exponentialpiston models (Maloszewski and Zuber, 1982;McGuire and McDonnell, 2006;McDonnell et al., 2010). For the WLT stream and river waters, Morgenstern and Taylor (2009) used LPM, which consists of two parallel exponential-piston flow models, to obtain mean residence times (MRTs). The tritium time series data for the Whareroa Streams, Kuratau and Omori rivers allowed the unique age interpretation with MRTs of 24, 12 and 18 yr, respectively. Waihaha and Whanganui rivers have tritium time series data that still have a degree of ambiguity in age interpretation, with MRT between 4 and 17 yr. Groundwater age distributions with MRTs for the five river catchments can also be obtained with the use of the tritium-calibrated MT3DMS model. Simulating groundwater flow and tracer transport (e.g. 3 H) with distributed parameter models (DPMs) allows the inclusion of a high level of complexity into the flow model unavailable in LPMs, e.g. spatially distributed groundwater recharge, surface water networks, and aquifer hydrogeology. Once calibrated to tritium data and spatial variability input histories of chemicals (e.g. nitrate) included, such models allow us to obtain realistic nitrate concentration at any stream reach. From the tritium-calibrated MT3DMS model, the groundwater age concentration was assigned as a single mobile species with negative zero-order decay reaction that allows aging of groundwater by a day per day (Zheng, 2009). Zero age was assigned to the top of the model (as recharge concentration), and the transient MT3DMS model was run for 200 yr to reach steady-state conditions of groundwater age concentrations. The simulated groundwater age concentrations were determined in each drain cell that is included in the surface water network of the river catchment, and the groundwater age distribution was constructed for the stream outlets of the five catchments. The simulated groundwater age distributions for the Waihaha, Whanganui, Whareroa, Kuratau and Omori surface water catchments are shown in Fig. 5 with MRTs of 9.5, 6.1, 15.1, 7.1 and 3.3 yr, respectively. In Fig. 5, these MRTs correspond to 55.1 %, 65.1 %, 49.7 %, 62.4 % and 56.4 % of the cumulative frequency distributions (CFDs) for the five rivers. This implies that the simulated CFDs, except for that at Whareroa Stream, are very similar to the exponential CFD shape and MRT of the exponential CFD. The MRT is equal to 63 % of the CFD for the exponential shape (Haitjema, 1995;Abrams, 2012). For Waihaha River and Whanganui Stream, the MRTs simulated with DPM are within the wide ambiguous MRT ranges estimated with LPM. The MRTs for Whareroa Stream, Kuratau and Omori rivers are unique age solutions from the LPM, and the estimated MRTs via the DPM are much younger than those estimated from LPM. The MRT discrepancies between LPM and DPM for the Whareroa Stream, Kuratau and Omori rivers may be due to the alternative conceptual representation of WLTC in the DPM, e.g. spatially distributed groundwater recharge, surface water networks, and aquifer hydrogeology. For example, the Waihaha and Whanganui catchments have relatively simple hydrogeology compared to Whareroa Stream, Kuratau and Omori river catchments (Figs. 1 and 2). However, this investigation of alternative conceptual models for the WLTC is beyond the scope of this study.
To understand these matches and discrepancies in MRTs, the simulated groundwater ages in the mobile domain are shown in Fig. 6 for several model cross-sections that correspond to the investigated catchments. These show complex surface-groundwater flow interactions influenced by aquifer hydrogeology. In cross-section 129 for the Waihaha catchment, the youngest water (< 5 yr) occurs near the water table and groundwater becomes older with aquifer depth (up to 50 yr at the bottom of the Whakamaru unit). Old groundwater is likely to upwell from the bottom of the aquifer and discharges into the stream network, which are referred to as strong sinks (Abrams, 2012). This results in a mixture of old and young groundwaters in the river water at several places in Fig. 5 and can be conceptualized as the exponential flow pattern (Maloszewski and Zuber, 1982). In Fig. 5, CFDs for Waihaha, Whanganui, Kuratau and Omori are similar to that of the exponential model in spite of the hydrogeologic settings of these catchments being much more complex (see Figs. 1 and 2). In cross-section 213, the effect of two spatially distributed hydrogeologic units with different aquifer properties is demonstrated for the Whanganui catchment. The groundwater age for the Pakaumanu unit (on the left) is very young (< 5 yr) and for the Whakamaru unit (on the right) is up to 20 yr old (Fig. 6). A similar "groundwater aging" effect of the Whakamaru unit (W5a) can be observed in crosssections 299, 330 and 386, where the simulations result in older ages in the areas of the Whakamaru unit compared to younger ages simulated in the surrounding units (see Figs. 2 and 6). A groundwater divide is also observed to the left of the vertical unit boundary, which makes the river behave as a weak sink that captures only part of the old groundwater arriving from the divide, while the rest passes underneath the river towards the lake (Abrams, 2012). A similar pattern is observed in cross-section 386 that shows the Whakamaru unit in the Kuratau and Omori river catchments. The age of water at the bottom of the aquifer below the Omori catchment is much younger than that in the Kuratau catchment. The groundwater divide occurs near the catchments' boundary and separates the two catchments hydraulically. Note that these groundwater divides do not always coincide with the boundaries of the surface water watersheds.
In cross-sections 299 and 330 for the Whareroa and Kuratau catchments, complex exponential and piston flow patterns are present as a result of surface-groundwater interactions and the presence of the more conductive Oruanui unit above the less conductive Whakamaru unit (W5b). For the Oruanui unit, the groundwater ages show an exponential flow pattern with young groundwater ages near the water table and ages increasing with aquifer depth. For the underlying Whakamaru unit, a piston flow pattern occurs due to groundwater flow from the groundwater divide at the Pakaumanu unit toward Lake Taupo and is shown by vertical groundwater age contour intervals. The piston flow continues underneath the Oruanui unit, and the old groundwater upwells into rivers at the surface; see the middle of cross-sections 299 and 330. In cross-section 330 for the Kuratau catchment, groundwater with young ages underflows the river from the surrounding groundwater watersheds in the Whakamaru and Pakaumanu units. This results in a mixture of young and old groundwater underneath the river and indicates that rivers and streams may receive old groundwater depending on the hydrogeologic conditions in the aquifer.

Concluding remarks
We have implemented a second calibration approach by Zuber (2011) of a transient MT3DMS transport model to measured tritium concentrations in river and stream waters. Steady-state MODFLOW groundwater flow model and transient MT3DMS transport model were developed for the western Lake Taupo catchment (WLTC) using an 80 m uniform model grid to include small surface water features. In our study, we are using off-the-shelf package US Geological Survey finite-difference code MODFLOW-2000 (Mc-Donald andHarbaugh, 1988;Harbaugh et al., 2000) to model groundwater levels and baseflows and finite-difference code MT3DMS v. 5.2 to model tritium concentrations in ground-water and baseflows. The MT3DMS model was developed by Zhang (1990) on US Corps of Engineers funding to simulate advection, dispersion, and chemical reactions of contaminants in groundwater systems. The results of the groundwater flow model gave groundwater elevations ranging from 357 m above mean sea level at the Lake Taupo lakeshore to over 1000 m in the northern part of the model domain. For our study area, we adopted a simple water balance procedure of groundwater recharge estimation for steady-state regional models without sufficient time series field data (Sanford, 2002). This steady-state assumption allowed us to investigate groundwater contributions to stream baseflows at the average climate conditions in the WLTC.
Tritium measurements for the Waihaha, Whanganui, Whareroa, Kuratau and Omori rivers of the WLTC were used to calibrate the tritium transport model while maintaining the groundwater flow model (MODFLOW) calibration to groundwater elevation and stream baseflow data. The Kuratau catchment had tritium results for 1960-1970 in addition to the 2000-2009 measurements. Tritium data from around the bomb peak are very sensitive to groundwater mixing processes and the mean residence time of the groundwater, and therefore allowed us to fine-tune the transport model. The transport model calibration to tritium measurements in streams was conducted by adjusting one parameter at a time with simple calibration techniques, due to Visual MOD-FLOW software limitations (SWS, 2010), and was challenging due to different and complex surface-groundwater interactions in the five river catchments of the WLTC. For groundwater flow and transport model calibration, sophisticated calibration techniques such as the parameter estimation software PEST (Doherty, 2010) could be used to provide insights into uncertainties of field data and model structure. However, time and expertise are required to set up a PEST run, and PEST simulations are expected to take considerable time, especially for the large times series datasets in transient models. In addition, fully coupled groundwater flow and tritium transport models may be required to simulate changes in tritium concentrations in large rivers due to river losses to groundwater. Currently, only a few groundwater flow and transport models have this capability.
Tritium concentrations and groundwater ages were simulated with single-and dual-domain models. The simulated tritium concentrations given by the dual-domain model matched the measured data best in both the 1960-1970 and 2000-2008 periods, while those for the single-domain model matched 1960-1970 measurements but were too low for 2000-2009 measurements. In the single-domain model, tritium concentrations occur in flowing groundwater, whereas stagnant groundwater is present with tritium concentrations that are exchanged with the mobile phase using the mass transfer coefficient in the MT3DMS model in the dualdomain model. High and low values of the mass transfer rate in the dual-domain model make it equivalent to the singledomain model with total porosity and the single-domain model with the porosity of the mobile domain, respectively (Zheng and Wang, 1999). Overall, the tritium calibration approach is applicable to both the single-and the dual-domain MT3DMS settings.
The tritium-calibrated MT3DMS model was applied to simulate the groundwater ages. Groundwater transit time distributions were constructed with MT3DMS for the Waihaha, Whanganui, Whareroa, Kuratau and Omori river catchments with MRTs of 9.5, 6.1, 15.1, 7.1 and 3.3 yr respectively. Morgenstern and Taylor (2009) used LPM, which consists of two parallel exponential-piston flow models, to determine ambiguous MRT between 4 and 17 yr for Waihaha River and Whanganui Stream and unique MRTs of 24, 12 and 18 yr for the Whareroa Stream, Kuratau and Omori rivers, respectively. For Waihaha River and Whanganui Stream, the MRTs simulated with DPM are within the wide ambiguous MRT ranges estimated with LPM. MRTs for Whareroa Stream, Kuratau and Omori rivers are unique age solutions from the LPM, and the estimated MRTs via the DPM are much younger than those estimated from LPM. The MRT discrepancies between LPM and DPM for the Whareroa Stream, Kuratau and Omori rivers may be due to the alternative conceptual DPM of WLTC, e.g. spatially distributed groundwater recharge, surface water networks, and aquifer hydrogeology, and was beyond the scope of this study.
Complex surface-groundwater flow interactions in the WLTC aquifers are influenced by aquifer hydrogeology and indicated by the groundwater ages, levels and velocities. For the Waihaha catchment, groundwater ages indicate exponential flow patterns with young water (< 5 yr) near the water table and water becoming older with depth (up to 50 yr at the aquifer bottom). The groundwater age at the aquifer bottom below the Omori River water catchments is much younger than that in the Kuratau catchment. This can be attributed to the groundwater divide that separates the two catchments hydraulically. For the Whanganui catchment, the groundwater age for the Pakaumanu unit (on the left) is very young (< 5 yr) and for the Whakamaru unit (on the right) is up to 20 yr old. A groundwater divide is also observed to the left of the vertical unit boundary and makes the river behave as a weak sink capturing only part of the old groundwater arriving from the divide while the rest passes underneath the river towards the lake. For the Whareroa and Kuratau catchments, complex exponential and piston flow patterns are present as the result of surface-groundwater interactions and the presence of the more conductive Oruanui unit above the Whakamaru unit. This indicates that rivers and streams may receive old groundwater depending on the hydrogeologic conditions in the aquifers. Therefore, the calibrated tritium model and groundwater transit time distributions illuminate complex surface-groundwater interactions and are potentially important for land use pollution management in the WLTC.
This study is the first attempt to use tritium measurements in streams and rivers for groundwater flow and transport model calibration and illustrates the benefits of tri-tium data. In future studies, the long-term effect of seasonal climate variability on groundwater contributions to stream baseflows can be investigated with the transient groundwater flow model. A considerable amount of field data would be required to construct and calibrate a transient groundwater flow model on the regional scale for the past and current climatic variability. A transient groundwater flow model would be the next step for this modelling exercise, but lack of field data did not allow us to construct a transient model at present. Transient groundwater flow and tritium transport models are currently being developed in different study areas in NZ.
The overall purpose of long-term research is to match measured tracer concentration with simulated tracer concentrations using DPMs and LPMs and to compare resulting groundwater ages in order to understand the effect of hydrogeology and surface water hydrology. The simple LPMs provide a quick and robust estimate of MRTs and are independent of surface water catchment scale. On the other hand, the numerical DPMs allow us to include a high level of hydrogeological complexity unavailable in LPMs. For example, the tritium-calibrated transport MT3DMS model for the WLTC includes surface water network, zones of groundwater recharge and several hydrogeological units. Crosscomparison between LPMs and DPMs allows us to better understand connections between hydrogeology and surface water hydrology, and to investigate alternative conceptual LPMs and DPMs of the WLTC to match their MRTs and groundwater age distributions.