Effectiveness of a Natural Headwater Wetland for Reducing Agricultural Nitrogen Loads

: Natural wetlands can play a key role in controlling non-point source pollution, but quantifying their capacity to reduce contaminant loads is often challenging due to diffuse and variable inflows. The nitrogen removal performance of a small natural headwater wetland in a pastoral agricultural catchment in Waikato, New Zealand was assessed over a two-year period (2011–2013). Flow and water quality samples were collected at the wetland upper and lower locations, and piezometers sampled inside and outside the wetland. A simple dynamic model operating on an hourly time step was used to assess wetland removal performance for key N species. Hourly measurements of inflow, outflow, rainfall and Penman-Monteith evapotranspiration estimates were used to calculate dynamic water balance for the wetland. A dynamic N mass balance was calculated for each N component by coupling influent concentrations to the dynamic water balance and applying a first order areal removal coefficient (k 20 ) adjusted to the ambient temperature. Flow and water quality monitoring showed that wetland was mainly groundwater fed. The concentrations of oxidised nitrogen (NO x -N, Total Organic Nitrogen (TON) and Total-N (TN) were lower at the outlet of the wetland regardless of flow conditions or seasonality, even during winter storms. The model estimation showed that the wetland could reduce net NO x -N, NH 4 -N, TON and TN loads by 76%, 73%, 26% and 57%, respectively. measured using a ﬁeld 15 N tracer addition approach.


Introduction
Nitrogen loads from diffuse water pollution (DWP) are a major water quality problem in many countries [1]. Compared to point source pollution, DWP is more complex and difficult to control due to its numerous and dispersed sources, and the difficulties in tracing its pathways [2]. Wetlands have been demonstrated to be an effective means to attenuate nitrogen derived from DWP [3,4] by plant and periphyton uptake and microbial denitrification [5,6]. In some wetlands, studies have shown that denitrification is the dominant nitrate removal mechanism [4,7,8]. However, microbial activity, which controls denitrification rates, can be reduced at low temperatures [9], pH, and carbon availability [6,10]. In line with these seasonal differences in wetlands, nitrogen removal performance is widely reported with reduced rates measured at colder temperatures [4]. Bernal [11] examined N content, N accumulation rates and soil C:N ratios over time in two riverine wetlands in the U.S. Midwest and found that, besides denitrification, organic accumulation was also important in N removal. Zaman [12] found that denitrification only accounted for 6-7% of observed NO 3 removal in a short-term injection-resampling study, and suggested that plant uptake was the principal removal mechanism.
In New Zealand, about half of the land area is used in some form of pastoral production, ranging from intensively farmed lowland (often dairy) to extensively farmed hill country (usually sheep/beef/deer) [13]. Dairy cow numbers have risen from approximately 3 million in 1980 to approximately 6.5 million in 2015 [14] along with increased application of fertilisers and use of supplementary feed. For example, N fertiliser use increased from 50 Gg in 1989 to 329 Gg in 2010 [15]. At the same time, extensive drainage across many parts of New Zealand have converted large areas of natural wetlands into agricultural land [3] and small remnant wetland areas continue to be drained as farming practices intensify [16]. In New Zealand, it has been estimated that over 90% of the former wetland area has been lost within a century and a half [17] and the trend is continuing particularly for small wetlands in agricultural landscapes. Ongoing intensification has raised concerns about environmental sustainability and contamination of groundwater and surface water with nutrients particularly under intensive dairying [18,19]. Therefore, there is increasing need for effective management tools to reduce the nutrient losses to water bodies [13].
Natural headwater wetlands are a relatively common feature in the hilly parts of New Zealand. These wetlands often occur within the headwater areas of catchments and along the sides of streams [20]. Although they are individually small, they may represent a significant proportion of headwater catchments. The potential of these wetlands to attenuate upslope derived pollutants is well recognised [21]. However, they can also be a potential source of agricultural pollutants because of their direct connection to the stream network, and farmers see them as a suitable drinking water source for livestock [20]. Pastoral wetlands are often small (less than 5000 m 2 ) and therefore they are rarely identified in wetland inventories or managed, and the ecosystem services they provide are poorly understood. Compared to constructed wetlands, the contaminant removal processes in natural wetlands are generally more complicated to study because they have diffuse, spatially and temporally variable groundwater inflows and outflows that are difficult to access and measure. There have been few attempts to quantify the effectiveness of natural seepage wetlands [22], whereas nutrient removal of constructed wetlands with discrete inflows and outflows is comparatively well studied [3,[23][24][25][26].
Models are being increasingly used to help quantify and predict the efficacy of wetlands for removing excess nutrients. The biological, chemical and physical removal processes in wetlands vary in space and time [11]. Modelling enables improved understanding of the complex interplay of hydrology and biogeochemical processes taking place in wetlands and their variability in space and time at relatively low cost.
This study aims to estimate the performance of a natural headwater wetland in removing nitrogen loads originated from upland grazed dairy pasture. It was hypothesised that wetland N removal efficiency would vary depending on inflow loads and seasonal temperatures, and that the coincidence of high flows and low temperatures in winter would reduce removal efficiencies. The performance of the wetland was estimated by assessing nitrogen in-loads and out-loads by flow and water quality monitoring and modelling.

Study Site
The study wetland is located on a dairy farm near Kiwitahi in the headwaters of the Toenepi catchment (15.8 km 2 ) in the eastern Waikato region, New Zealand ( Figure 1). The Toenepi catchment is intensively farmed and 75% of the catchment area is under dairy production with a stocking rate of 3 cows ha −1 [27]. The mean annual rainfall of the area is 1377 mm. The upper Toenepi catchment is hilly with~80% of the area classified as either rolling or steep (>10% gradient) [28]. The wetland catchment is dominated by Morrinsville clay soil (NZ Soil Classification: Orthic Granular). Rotational grazing of a single herd of ~220 (Holstein Friesian) cattle is practiced on the dairy farm, which is divided into 33 individual paddocks or fields (fenced pasture area for grazing) of between 1.0 and 3.1 ha. The study wetland is located within a small (~1.9 ha) fenced paddock, which is grazed for ~1 day every 40 days during winter and summer and ~1 day every 20 days during spring and autumn. As most of the time the wetland does not provide ready access to surface water, drinking water is available to the herd from a water trough (groundwater bore source) within the paddock. The paddock containing the study wetland (with exception of the wetland itself) is steep (mostly exceeding 20° slope). The extent and impacts of grazing events on water quality at the outlet of the wetland over the period of the present study have been described by Hughes [29].
The wetland has an area of ~0.15 ha (2.8% of surface catchment) and an average slope of 3.5°. The permanently saturated wetland soil is composed largely of a mix of organic material and the clay-based soils eroded from the surrounding hillslopes. Sediment probe measurements indicated that within 1 m of the edge the saturated layer was generally between 0.5 and 1 m thick. Depths increased with distance from the edge, and were generally between 1 and 2 m in the centre of the wetland. The wetland vegetation is dominated by glaucous sweet grass (Glyceria declinata), a perennial aquatic grass widely naturalised in New Zealand.

Site Monitoring
The study site was monitored for a two-year period between October 2011 and September 2013. Flow was measured hourly at two 45° v-notch weirs with stage height measured by a Unidata Hydrologger water level recorder (1 mm resolution; Unidata Pty, O'Connor, WA, Australia) and converted to flow using a theoretical rating equation. The Upper Weir (UW) was located at the head of wetland, downstream from where there was significant ground water seepage ( Figure 1). The catchment area above the UW is ~2.9 ha. The Lower Weir (LW) was located within a constricted part of the lower wetland with a total catchment area of ~5.2 ha (Figure 1). During the period between November 2012 and May 2013, the study site experienced exceptionally dry conditions [30]. Consequently, no flow was recorded at the UW from early January 2013 through to May 2013. Despite this, the wetland remained wet and boggy with flow at the LW.
An ISCO 3700 automatic water sampler (Teledyne Isco, Lincoln, NE, USA) was programmed to collect water samples behind each weir using a stage-based trigger. In addition, low flow grab samples were collected during site visits approximately every six weeks. Once collected, samples Rotational grazing of a single herd of~220 (Holstein Friesian) cattle is practiced on the dairy farm, which is divided into 33 individual paddocks or fields (fenced pasture area for grazing) of between 1.0 and 3.1 ha. The study wetland is located within a small (~1.9 ha) fenced paddock, which is grazed for~1 day every 40 days during winter and summer and~1 day every 20 days during spring and autumn. As most of the time the wetland does not provide ready access to surface water, drinking water is available to the herd from a water trough (groundwater bore source) within the paddock. The paddock containing the study wetland (with exception of the wetland itself) is steep (mostly exceeding 20 • slope). The extent and impacts of grazing events on water quality at the outlet of the wetland over the period of the present study have been described by Hughes [29].
The wetland has an area of~0.15 ha (2.8% of surface catchment) and an average slope of 3.5 • . The permanently saturated wetland soil is composed largely of a mix of organic material and the clay-based soils eroded from the surrounding hillslopes. Sediment probe measurements indicated that within 1 m of the edge the saturated layer was generally between 0.5 and 1 m thick. Depths increased with distance from the edge, and were generally between 1 and 2 m in the centre of the wetland. The wetland vegetation is dominated by glaucous sweet grass (Glyceria declinata), a perennial aquatic grass widely naturalised in New Zealand.

Site Monitoring
The study site was monitored for a two-year period between October 2011 and September 2013. Flow was measured hourly at two 45 • v-notch weirs with stage height measured by a Unidata Hydrologger water level recorder (1 mm resolution; Unidata Pty, O'Connor, WA, Australia) and converted to flow using a theoretical rating equation. The Upper Weir (UW) was located at the head of wetland, downstream from where there was significant ground water seepage ( Figure 1). The catchment area above the UW is~2.9 ha. The Lower Weir (LW) was located within a constricted part of the lower wetland with a total catchment area of~5.2 ha (Figure 1). During the period between November 2012 and May 2013, the study site experienced exceptionally dry conditions [30]. Consequently, no flow was recorded at the UW from early January 2013 through to May 2013. Despite this, the wetland remained wet and boggy with flow at the LW.
An ISCO 3700 automatic water sampler (Teledyne Isco, Lincoln, NE, USA) was programmed to collect water samples behind each weir using a stage-based trigger. In addition, low flow grab samples were collected during site visits approximately every six weeks. Once collected, samples were immediately placed in an insulated storage bin containing an ice slurry. Samples were delivered to the NIWA-Hamilton Water Quality Laboratory on the day of collection for grab samples and within 24 h for samples collected from the automatic sampler.
Piezometers were installed at 13 locations, both within (nine piezometers) and adjacent to the wetland (four piezometers) ( Figure 1). The within-wetland piezometers were installed to below the depth of the deposited wetland material (i.e., into the weathered underlying bedrock material) and were slotted to collect water from depths of 1-2.5 m. During site visits for grab samples, piezometer water level was recorded and water samples were collected. Two overland flow (OLF) samplers were also monitored at the site, one in the gully immediately upstream of UW (OLF1) and one at the base of a swale on the northern margin of the wetland (OLF2).
Rainfall was also measured at the lower weir at 10 min intervals by an Ota tipping bucket rain gauge (Ota Keiki Seisakusho, Tokyo, Japan). Potential evapotranspiration was calculated using the FAO Penman-Monteith equation [31].

Laboratory Analysis
Collected water samples were analysed for, oxidised N (nitrite-and nitrate-N, here-after referred to as NO x -N), ammonium-N (NH 4 -N), Total Nitrogen (TN). Total Organic Nitrogen (TON) was calculated by subtraction (TN minus (NO x -N plus NH 4 -N)). A Lachat flow injection analyser (Hach, Loveland, CO, USA) was used for NO x -N, NH 4 -N (detection limit 1 mg/m 3 ), and for TN (detection limit 10 mg/m 3 ). All water samples were filtered after subsampling for TN (as well as total suspended solids and total phosphorus not reported here) with a Millipore ® syringe and filter holder containing a GF/C glass fibre pre-filter (47 mm diam., 1.2 µm pore size), and a Sartorius ® cellulose acetate membrane filter (47 mm diam., 0.45 µm pore size).
To determine whether pollutant concentrations varied by season, the Kruskal-Wallis statistical test was used to test for differences in median seasonal pollutant concentrations. The Kruskal-Wallis test can be considered to be a non-parametric version of the ANOVA test. The Mann-Whitney U test was used for comparing upper and lower weir median pollutant concentrations. A significance level of p < 0.05 was adopted in all tests. The statistical software package Statistica 12 (StatSoft; Tulsa, OK, USA) [32] was used to perform these tests.

Modelling Nitrogen Removal
We used N measurements and flow measurements as inputs to a simple dynamic model [32,33] to explore potential wetland N removal performance during 2012 when continuous flow data were available. The model was set up in Microsoft Excel (365 Pro Plus ver. 1708, Redmond, WA, USA) ( Figure 2) and uses Euler integration (Equation (1)). To simulate the internal hydraulic dynamics of the wetland, we used a five tanks-in-series approach (Figures 1 and 2) with areas from 56 m 2 to 637 m 2 ( Table 1). Tank 1 was situated above UW and Tank 5 terminated at LW. Each tank had at least one piezometer located near its outlet. Hourly measurements of flow at UW and LW, rainfall, and Penman evapotranspiration estimates were used to calculate a dynamic water balance for the wetland.  The following differential equation (Kadlec, 2012) was used for the dynamic water balance: where is time (h), + 1 refers to Tanki+1, refers to Tanki, is the volume of water in the tank (m 3 ), is the flow (m 3 /h), is the surface area for the tank (m 2 ), is the precipitation (m/h), is the evapotranspiration (m/h), is the surface runoff inflow (m 3 /h), and is the net groundwater inflow (m 3 /h). When weir is not present, then Qi is calculated from water balance based on inputs/outputs from neighboring tanks, precipitation, surface runoff and evapotranspiration calculated based on each tank's area.
Storm events above a certain threshold (rainfall intensity ≥ 1 mm/h, flow at weirs increasing and storm length ≥ 3 h) were assumed to always result in surface runoff. Altogether, 37 such storms were defined over the period of this study (12 months, January 2012 to December 2012), the shortest being 3 h and the longest 27 h. SR into each tank was calculated as P-ET during the storm event and multiplied by the catchment area of the tank. The subsurface infiltration rate of groundwater into or out of the wetland was derived from the water balance calculated at the outlet of the wetland, assuming that the wetland is a closed system. Changes in wetland depth and hence storage volume were simulated for a 45° V-notch weir set 35 cm above the base of the wetland using a modification of the Francis weir formula [34].
A dynamic N mass balance was calculated by coupling influent concentrations to the dynamic water balance and applying a first order areal removal rate coefficient (k) (Equation (2)) [35].
The mass balance for a tank is given by the differential equation: where C is concentration (mg N/m 3 ) and k is the removal rate coefficient (m/h). The removal rate coefficient (k) is well known to be temperature sensitive [35]. The modified Arrhenius equation was used for adjusting k to the ambient temperature with temperature factors derived from field-scale wetland studies (Table 2) [35,36]. The annual water temperature regime was The following differential equation (Kadlec, 2012) was used for the dynamic water balance: where t is time (h), i + 1 refers to Tank i+1 , i refers to Tank i , V is the volume of water in the tank (m 3 ), Q is the flow (m 3 /h), A is the surface area for the tank (m 2 ), P is the precipitation (m/h), ET is the evapotranspiration (m/h), SR is the surface runoff inflow (m 3 /h), and GW is the net groundwater inflow (m 3 /h). When weir is not present, then Q i is calculated from water balance based on inputs/outputs from neighboring tanks, precipitation, surface runoff and evapotranspiration calculated based on each tank's area. Storm events above a certain threshold (rainfall intensity ≥ 1 mm/h, flow at weirs increasing and storm length ≥ 3 h) were assumed to always result in surface runoff. Altogether, 37 such storms were defined over the period of this study (12 months, January 2012 to December 2012), the shortest being 3 h and the longest 27 h. SR into each tank was calculated as P-ET during the storm event and multiplied by the catchment area of the tank. The subsurface infiltration rate of groundwater into or out of the wetland was derived from the water balance calculated at the outlet of the wetland, assuming that the wetland is a closed system. Changes in wetland depth and hence storage volume were simulated for a 45 • V-notch weir set 35 cm above the base of the wetland using a modification of the Francis weir formula [34].
A dynamic N mass balance was calculated by coupling influent concentrations to the dynamic water balance and applying a first order areal removal rate coefficient (k) (Equation (2)) [35]. The mass balance for a tank is given by the differential equation: where C is concentration (mg N/m 3 ) and k is the removal rate coefficient (m/h). The removal rate coefficient (k) is well known to be temperature sensitive [35]. The modified Arrhenius equation was used for adjusting k to the ambient temperature with temperature factors derived from field-scale wetland studies (Table 2) [35,36]. The annual water temperature regime was simulated using a sinusoidal function calibrated to air temperatures measured at the site, as described by Kadlec [33]. We adjusted within the range of k 20 coefficients documented for surface flow wetlands ( Table 2). For modelling net TON and TN removal, an irreducible background concentration C* was used ( Table 2; Equation (2)). Surface runoff was assigned higher NH 4 -N, NO x -N, TON and TN concentrations than UW, based on limited measurements of surface runoff around the wetland ( Table 3). The median value recorded for each tank was used as constant input value for groundwater NH 4 -N, NO x -N, TON and TN (Table 3). Atmospheric deposition of dissolved inorganic forms of N in New Zealand is considerably lower than commonly found in continental Northern Hemisphere regions [33,37], ranging from 1-5 g N ha −1 year −1 [38]. Therefore, a very low constant input value (10 mg/m 3 ) for all forms of nitrogen from rainfall was used. To improve the stability of the model, we decreased the time-step during the storms to 10 min. Wetland performance was calculated as: and areal mass load removal rate was calculated: Areal mass load removal rate (mg/m 2 /day) = (Inload − Outload)/wetland area/inflow days. (4)

Hydrology of the Wetland
Annual precipitation of 1056 mm was recorded for the wetland catchment in 2012. About 40% of the rainfall occurred during winter, 30% in spring, 20% in summer and 10% in autumn ( Figure 3). During 2012, flow occurred on all days at the LW (lower weir) but only on 122 days at the UW (upper weir). During that 122-day period, there was 1836 m 3 inflow to the wetland as measured at the UW and 12,439.6 m 3 outflow from the wetland as measured at the LW. Flow measured at the upper weir accounted for less than 15% of the flow exiting the lower weir. This is despite the area above the upper weir contributing over half the catchment area and including a highly convergent gully immediately upstream from the wetland. This suggests that a large proportion (85%) of the flow exiting the wetland is derived from subsurface inflows entering between the upper and lower weirs.

Water Quality of the Wetland
Altogether, 56 water quality samples from the upper weir and 95 samples from the lower weir were analysed (Table 3, Figure 4). Baseflow samples were collected throughout the study period and flow events were sampled over a range of both summer and winter conditions. Water quality data for storms has been grouped together into summer/autumn (December-May) and winter/spring (June-November).
When comparing upper and lower weir, all the measured variables showed lower (NOx-N and TN statistically significantly, p < 0.05) concentrations at the lower weir, regardless of flow conditions or seasonality (Figure 4, Table 2). During the winter/spring storms, all measured variables had significantly higher values (p < 0.05) at the upper weir (Tables 3 and 4).
In terms of flow conditions, the concentration of all variables was highest during storm events (Figure 4). At the upper weir, all stormflow variables were highest during summer, which could be due to more frequent stock access to the wetland paddock during this period [29]. However, at the lower weir, all variables, except NH4-N, were highest for winter storms.

Water Quality of the Wetland
Altogether, 56 water quality samples from the upper weir and 95 samples from the lower weir were analysed (Table 3, Figure 4). Baseflow samples were collected throughout the study period and flow events were sampled over a range of both summer and winter conditions. Water quality data for storms has been grouped together into summer/autumn (December-May) and winter/spring (June-November).
When comparing upper and lower weir, all the measured variables showed lower (NO x -N and TN statistically significantly, p < 0.05) concentrations at the lower weir, regardless of flow conditions or seasonality (Figure 4, Table 2). During the winter/spring storms, all measured variables had significantly higher values (p < 0.05) at the upper weir (Tables 3 and 4).
In terms of flow conditions, the concentration of all variables was highest during storm events (Figure 4). At the upper weir, all stormflow variables were highest during summer, which could be due to more frequent stock access to the wetland paddock during this period [29]. However, at the lower weir, all variables, except NH 4 -N, were highest for winter storms.
TN statistically significantly, p < 0.05) concentrations at the lower weir, regardless of flow conditions or seasonality (Figure 4, Table 2). During the winter/spring storms, all measured variables had significantly higher values (p < 0.05) at the upper weir (Tables 3 and 4).
In terms of flow conditions, the concentration of all variables was highest during storm events (Figure 4). At the upper weir, all stormflow variables were highest during summer, which could be due to more frequent stock access to the wetland paddock during this period [29]. However, at the lower weir, all variables, except NH4-N, were highest for winter storms.   Fifty-nine groundwater samples were analysed. Fifteen of those were collected adjacent to the wetland area and 44 within the wetland. Only six discreet overland flow samples were collected from two sites, but, during storm-flows, OLF represented most of the flow at the upper weir. This confirmed that OLF entering the wetland had the most elevated levels of all measured pollutants ( Figure 5). Groundwater around the wetland always had relatively low N concentrations. In case of NOx-N, it was significantly lower (p < 0.05) than overland flow and piezometers sampled within the wetland. Despite the similarity of the nitrate concentrations from around the wetland, the NOx-N concentrations from the sample sites within the wetland showed remarkable spatial variation. Above the UW, the measured NOx-N concentrations were consistently higher than 3000 mg/m 3 . Higher NOx-N concentrations were also detected in piezometer NP3, which was sited where another gully entered the wetland from north downstream of the UW. NH4-N and TON concentrations were significantly lower in the wetland groundwater than in OLF.  Fifty-nine groundwater samples were analysed. Fifteen of those were collected adjacent to the wetland area and 44 within the wetland. Only six discreet overland flow samples were collected from two sites, but, during storm-flows, OLF represented most of the flow at the upper weir. This confirmed that OLF entering the wetland had the most elevated levels of all measured pollutants ( Figure 5). Groundwater around the wetland always had relatively low N concentrations. In case of NO x -N, it was significantly lower (p < 0.05) than overland flow and piezometers sampled within the wetland. Despite the similarity of the nitrate concentrations from around the wetland, the NO x -N concentrations from the sample sites within the wetland showed remarkable spatial variation. Above the UW, the measured NO x -N concentrations were consistently higher than 3000 mg/m 3 . Higher NO x -N concentrations were also detected in piezometer NP3, which was sited where another gully entered the wetland from

Modelling Wetland N Species Removal Efficiency
NOx-N was the dominant form of N (60.7% of in load) entering the wetland ( Figure 6) mainly via groundwater seepage (90.7%). TON also comprised a significant proportion of the N loading (29.9% of in loads) with NH4-N only 9.4% of the load. TON and NH4-N were also mainly entering via groundwater (76.8% and 95.1%, respectively). However, a significant portion of TON was also entering the wetland via surface runoff (23%). TON was the dominant form leaving the wetland (56.1% of out load) followed by NOx-N (37.5% of out load) and NH4-N (6.4%). Best fit k20 values based on least squares varied for nitrogen forms. NOx-N and TN showed best fit at the upper boundary of commonly recorded k20 rates (168 and 40 m year −1 , respectively; Table 2) and TON at lowest levels (5 m year −1 ). Calculated wetland removal efficiency (Equation (3)) varied considerably for different nitrogen forms (Table 5). Modelling results show substantial reductions in NH4-N and NOx-N loads after passage through the wetland, except during the winter period when load reductions were more muted (Figure 7a,b; Table 5). Although the overall NH4-N and NOx-N mass removal efficiency for the wetland (Equation (3)) was 72.9% and 75.5%, respectively, the   (Figure 6) mainly via groundwater seepage (90.7%). TON also comprised a significant proportion of the N loading (29.9% of in loads) with NH 4 -N only 9.4% of the load. TON and NH 4 -N were also mainly entering via groundwater (76.8% and 95.1%, respectively). However, a significant portion of TON was also entering the wetland via surface runoff (23%). TON was the dominant form leaving the wetland (56.1% of out load) followed by NO x -N (37.5% of out load) and NH 4 -N (6.4%).

Modelling Wetland N Species Removal Efficiency
NOx-N was the dominant form of N (60.7% of in load) entering the wetland ( Figure 6) mainly via groundwater seepage (90.7%). TON also comprised a significant proportion of the N loading (29.9% of in loads) with NH4-N only 9.4% of the load. TON and NH4-N were also mainly entering via groundwater (76.8% and 95.1%, respectively). However, a significant portion of TON was also entering the wetland via surface runoff (23%). TON was the dominant form leaving the wetland (56.1% of out load) followed by NOx-N (37.5% of out load) and NH4-N (6.4%). Best fit k20 values based on least squares varied for nitrogen forms. NOx-N and TN showed best fit at the upper boundary of commonly recorded k20 rates (168 and 40 m year −1 , respectively; Table 2) and TON at lowest levels (5 m year −1 ). Calculated wetland removal efficiency (Equation (3)) varied considerably for different nitrogen forms (Table 5). Modelling results show substantial reductions in NH4-N and NOx-N loads after passage through the wetland, except during the winter period when load reductions were more muted (Figure 7a,b; Table 5). Although the overall NH4-N and NOx-N mass removal efficiency for the wetland (Equation (3)) was 72.9% and 75.5%, respectively, the  Table 2) and TON at lowest levels (5 m year −1 ). Calculated wetland removal efficiency (Equation (3)) varied considerably for different nitrogen forms (Table 5). Modelling results show substantial reductions in NH 4 -N and NO x -N loads after passage through the wetland, except during the winter period when load reductions were more muted (Figure 7a,b; Table 5). Although the overall NH 4 -N and NO x -N mass removal efficiency for the wetland (Equation (3)) was 72.9% and 75.5%, respectively, the removal efficiency during the winter storms was only 43.5% for NH 4 -N and 67.3% for NO x -N. This is likely to result from the coincidence of high flows (resulting in reduced hydraulic residence times; HRT) and low water temperatures during winter (Figure 7d). However, despite reduced removal efficiency NO x -N and NH 4 -N mass removal was actually greater over the winter period because of the higher and more consistent loads received ( Table 5). The modelling showed only 26% of net TON reduction in the wetland, with highest rates during summer/autumn storm events. Annual reduction of TN was 57.2%, varying from 39.8% during baseflow to 75% during summer/autumn storms.
Water 2018, 10, x FOR PEER REVIEW 10 of 18 removal efficiency during the winter storms was only 43.5% for NH4-N and 67.3% for NOx-N. This is likely to result from the coincidence of high flows (resulting in reduced hydraulic residence times; HRT) and low water temperatures during winter (Figure 7d). However, despite reduced removal efficiency NOx-N and NH4-N mass removal was actually greater over the winter period because of the higher and more consistent loads received ( Table 5). The modelling showed only 26% of net TON reduction in the wetland, with highest rates during summer/autumn storm events. Annual reduction of TN was 57.2%, varying from 39.8% during baseflow to 75% during summer/autumn storms.

Discussion
Water quality monitoring showed that the concentrations of NH 4 -N, NO x -N, Total Dissolved Nitrogen and Total-N were consistently lower at the outlet of the wetland regardless of flow conditions or seasonality, and, during the winter storms, all pollutants had significantly lower values. This was a strong indication that the wetland is very efficient at attenuating the poor water quality entering from upstream. Although the water quality measurements at the upper and lower wetland weirs indicated highest concentrations for TON in the surface water, the N entering in groundwater was mainly in the form of NO x -N, with relatively lower TON concentrations. The dynamic hydrological sub-model calibrated with flow data from upper and lower weirs and local climatological data provided a rational approach to estimate relative inflows to the wetland from different pathways. It showed the wetland was almost 95% groundwater fed. Linking this information with median concentrations of N species measured in the different inflows to the wetland showed that NO x -N was the dominant form of N entering (estimated as~61% of load). The first order removal model estimated annual wetland net removal efficiency to be approximately 76% for NO x -N, 73% for NH 4 -N and 26% for TON. These forms of N are of course subject to transformations between different organic and inorganic forms via microbial processes (nitrification, denitrification, anammox and DNRA), assimilation by plants and microbes, adsorption to particulate phases (particularly NH 4 -N), and settling and storage of particulate TON in the wetland. We did not have sufficient data about relative changes in N species or balance between different processes to be able to usefully infer sequential processing rates (e.g., [39]) in these large and highly dynamic systems.
The apparent treatment efficiency (percentage removal) of the headwater wetland was high for NO x -N and TN compared with that of surface-flow constructed wetlands receiving tile drainage [3,33]. The wetland influent N loading rates in the present study were relatively low (TN~60, NO x -N 33 mg m −2 day −1 ) compared to rates in excess of 400 mg m −2 day −1 measured in the constructed wetlands receiving tile drainage [3]. However, the best fit k 20 reaction rate coefficient for NO x -N was also very high compared to those measured in field scale constructed wetlands [4,33].
As with many headwater pastoral headwater wetlands in New Zealand, our study wetland is likely to have formed since European settlement of New Zealand. Since the catchment was cleared of its indigenous forest land cover during the 19th and early 20th centuries, clay-rich soils have eroded from the steep slopes and accumulated in the bottom of the gully to form a wetland swale stabilised by the dense sward of amphibious glaucous sweetgrass (Glyceria declinata). Similar features have been described from south-east Australia [40]. Our study wetland was strongly was by groundwater. However, the wetland was also losing water occasionally to groundwater. We detected slight spatial declining trend along the wetland in nitrate-N and Total-N concentrations, which has been also detected by other studies. Xu et al. [9] also detected decline in Total-N and nitrate-N concentrations along the wetland but not in the concentrations of ammonium, which may be removed by transient sorption, nitrification, anammox and plant uptake, but also generated through mineralization of organic N, and dissimilatory nitrate reduction to Ammonium [4]. Ackerman et al. [41], Maxwell et al. [42], and Brauer et al. [43] detected significant pollutant reduction along the groundwater flow paths associated with seepage from a wetland. They illustrate that subsurface flux away from wetlands provides another mechanism for N removal. Shallow groundwater infiltrating though the base and sides of the wetland would pass through saturated layers of anoxic sedimentary and organic material and plant detritus with high denitrification potential [7,35,44]. High denitrification rates are common in wetlands and riparian areas that are rich with organic material. For example, Whitmire and Hamilton [45] using push-pull techniques found that nitrate added to groundwater and re-injected into wetland sediments disappeared rapidly relative to conservative tracers without any lag time, and was depleted to below detection limits (10-15 µg N L −1 ) within 5 to 20 h. In a 15 N denitrification study in a similar natural riparian wetland intercepting shallow groundwater, Burns and Nguyen [46] measured >90% removal of injected 15 NO 3 -N along a 100 cm flow path through the organic wetland soil, with essentially all of the 15 NO 3 -N-removed within 30 cm. Extensive 15 N tracer study performed on small headwater agricultural streams in the USA by Peterson [21] showed that NH 4 -N and NO x -N entering the streams was removed within tens to hundreds of meters and in one of the streams the denitrification rate was 0.045 µg N m −2 s −1 [47]. Denitrification consisted almost entirely of N 2 production (99%), with very little N 2 O production occurring.
To describe the temperature response of denitrification, we used Arrhenius equation, which is most commonly used in this case and where the temperature has a stronger effect on the denitrification at lower temperatures [4,47]. Several authors have reported that denitrification has a "breakpoint" temperature, below which denitrification rates decrease more rapidly [9,48,49]. The exact temperature of the breakpoint is not known, but, in most of the studies, it varies usually between 10 and 12 • C [50]. However, Holtan-Hartwig et al. [51] found significant decrease in denitrification rates below 0 degrees and reported major deviations from the regular Arrhenius response at near freezing temperatures. Our results showed lower removal efficiency between 8 and 10 degrees, which coincided with high flows from surface runoff, which also had high loads of pollutants. Some studies have showed that wetlands receiving very variable flows show reduced nitrate removal performance [3,52], whereas Hernandez and Mitsch [53] found in the opposite that denitrification rates in the high marsh zone were significantly higher under flood pulsing (778 ± 92 µg N m −2 h −1 ) than under steady flow (328 ± 63 µg N m −2 h −1 ), but, in the low marsh and edge zones, flood pulses did not affect denitrification. Kadlec [4] also argues that, if wetland receives pollutants in pulses and the water interval between flows is sufficiently long, then the wetland will contain very low nitrate concentrations when another pulse arrives due to inter-event treatment. In the present study, removal efficiency declined during winter period, which coincided with low temperatures. Plants have been found to improve the wetland removal efficiency during the cold season. Yates et al. [54] demonstrated that Carex aquatilis helped to enhance nitrogen removal at cold temperatures and high hydraulic load. In addition, extended HRT could improve TN removal at low temperatures [9].
High denitrification rates are likely to be stimulated by the Glyceria vegetation that dominated the wetland. Matheson et al. [8], investigating the fate of 15 N-NO 3 in wetland soil microcosms, found that, in the presence of Glyceria, denitrification was elevated compared to unplanted microcosms and was the principal mechanism of NO x -N removal (61-63%). They attributed high denitrification rates and low Dissimilatory Nitrate Reduction to Ammonium (DNRA) to differences in soil redox state induced by the presence of the plant. During the growth season, plant and periphyton uptake of dissolved inorganic forms of N is also likely to have occurred. Subsequent senescence and decay of this plant material is likely to have generated autochthonous particulate and dissolved organic N that was exported from the wetland.
Organic N export from agricultural land is often overlooked as a significant pathway of diffuse N loss from pasture lands [55,56], which, in some instances, can be of similar magnitude to NO x -N losses. The estimated load of TON exiting from our wetland was 6.5 kg, which equates to 1.3 kg ha −1 . This is lower than the median level reported by van Kessel [55] (4.0 kg N ha −1 year −1 ), and compared to dissolved organic N leaching rates measured in lowland pastoral soils elsewhere in New Zealand of 28-117 kg N ha −1 year −1 [57].
TON, the second largest component of TN entering the wetland, passed through the wetland with 26% net reduction. Wetlands characteristically accumulate a large pool of organic C and N, and experience with constructed wetlands suggests incoming organic nitrogen is reduced to a non-zero background concentration, created by residuals and wetland return fluxes [35]. A significant proportion of the TON load entered the wetland via surface runoff during storms, suggesting that much may have been associated with particulates [20]. Wilcock [36] previously studied the narrow (1-3 m) riparian-planted wetland swale that starts~65 m below the lower weir of the headwater wetland in the present study. Based on their first two years of data, it is likely that around 50-60% of the TON exiting the headwater wetland would have been particulate organic N, likely associated with fine particulate detritus from plant litter and sloughed biofilm flocs. The remainder of the TON exported is expected to be dissolved organic forms, a significant proportion of which is likely to be bioavailable [58,59]. Analysis of the mass loading and removal data of Wilcock [36] shows that around 80% of the DON was removed in the~350 m length of relatively aerobic wetland swale downstream of our lower weir, suggesting that much of it was not recalcitrant. There are several processes that can lead to a decrease in the DON concentration as it moves through the soil profile. Dissolved organic N is used as a substrate by soil microbes [55]. Wymore et al. [60] hypothesize that under a high N scenario microbes may be more energy limited (i.e., low C:N ratio) and DON may serve primarily as an energy source. However, they also found strong underlying seasonal patterns in the response to NO 3 addition, indicating that the role of DON can switch between serving as a nutrient source to an energy source. The other process that can lead to a decrease in the DON concentration is the uptake of DON by plants. The uptake of organic N has been found to be widespread across different ecosystems and to consist of an important source of N, in particular in ecosystems where microbial biomass is prone to large seasonal fluctuations and contributes to the release of labile organic N [55,61]. For modelling of TON removal, we used a relatively low background concentration (20 mg/m 3 ) based on concentrations measured during very low flow periods and in the wetland downstream [36]. The effect of background concentration on modelling results was evident as most of the time wetland received TON levels only slightly higher than the background concentration and therefore the first order removal was applied to a very small proportion of TON, which resulted also in low removal efficiency. During the storms, TON in loads was considerably higher than background concentration and the removal efficiency was also higher. There also appears to be little or no temperature dependence of organic nitrogen k-values [35]. This was also evident in this wetland as summer and winter TON concentrations for lower weir did not differ significantly (Figure 4). For modeling, we used a low temperature coefficient (1.01, Table 2), which resulted in very small seasonal variation in TON removal. NH 4 -N constituted only 9.4% of the inflow and the overall reduction was up to 73%. Most of the NH 4 -N likely originated from the mineralised urea in cattle urine patches. The amount of N under a urine patch can be equivalent to an application rate of 700 to 1200 kg N ha -1 [62], much higher than the demand of N for any agricultural crop [55]. NH 4 -N is commonly the preferred form of N for uptake by wetland plants, and, therefore, a proportion may have been assimilated by plants within the wetland during the growth season.
It is interesting to note that cattle grazing of the study wetland had only a very limited immediate impact on water quality [29]. According to Hughes et al. [29], a measurable increase in pollutants (attributable to cattle generated disturbance) occurred only once out of 18 high-intensity grazing days. This occurred when a cow became entrapped in close proximity to the downstream monitoring weir on a day when wetland flow was elevated.
The nitrogen (especially NO x -N and NH 4 -N) removal efficiency of our wetland was likely reduced by variability of inflow rates because the wetland received high nutrient level pulses during the winter storms, which coincided with lower temperatures [33]. This may have been exacerbated by water flowing across the wetland surface and around the wetland edges during high flow events [22,44]. Although the removed load during the storms was on an average 10 times higher than during baseflow, the removal efficiency (percentage removal) was lower. Nutrient removal by the wetland could be further improved by increasing residence time in the wetland through construction of a series of check dams and/or outflow control structures, although these would be prone to damage and erosion during storm-flows. The alternative of buffering inflow variability would appear to be practically difficult in the steep terrain surrounding the wetland, without changing from pasture to shrub or forest.

Conclusions
In comparison to constructed wetlands receiving defined surface flows (e.g., wastewaters or tile drainage), estimating the nutrient removal effectiveness of natural wetlands receiving a variable mix of seepage and overland flows is challenging. The measured concentrations of nitrate-N, Dissolved Inorganic Nitrogen and Total-N were always lower at the outlet of the wetland regardless of flow conditions or seasonality, even during winter storms. This indicated that the wetland is an effective removal mechanism for N. Use of the simple dynamic model calibrated to continuous flow monitoring, meteorological data and periodic water quality sampling enabled the loads from different pathways to be partitioned and load reductions estimated. Despite the highly variable flow of the wetland, it showed very high net overall reduction of NO x -N and NH 4 -N. Seepage of shallow nitrate-rich groundwater in through the organic-rich, anaerobic base of the wetland is likely to have stimulated the high denitrification rates modelled. Net removal of TON was less effective, likely due to decay processes in the wetland ecosystem, which regenerate N from plant detritus and biofilms, resulting in a natural background concentration. The results show that small seepage wetlands in the headwaters of New Zealand streams can be very effective at removing nitrogen loads. As their nutrient removal efficiency per unit land area is high, farmers might be more willing to maintain these natural wetlands undrained and fence them to prevent stock access. The study helps to raise awareness of the functions and values of these small pastoral seepage wetlands located in the headwater areas of the catchments.