Biogeosciences Modelling the effect of aggregates on N 2 O emission from denitrification in an agricultural peat soil

Nitrous oxide (N2O) emissions are highly variable in time, with high peak emissions lasting a few days to several weeks and low background emissions. This temporal variability is poorly understood which hampers the simulation of daily N2O emissions. In structured soils, like clay and peat, aggregates hamper the diffusion of oxygen, which leads to anaerobic microsites in the soil, favourable for denitrification. Diffusion of N2O out of the aggregates is also hampered, which leads to delayed emissions and increased reduction of N2O to N2. In this model simulation study we investigate the effect of aggregates in soils on the N 2O emissions. We present a parameterization to simulate the effects of aggregates on N 2O production by denitrification and on N2O reduction. The parameterization is based on the mobileimmobile model concept. It was implemented in a field-scale hydrological-biogeochemical model combination. We compared the simulated fluxes with observed fluxes from a fertilized and drained peat soil under grass. The results of this study show that aggregates strongly affect the simulated N2O emissions: peak emissions are lower, whereas the background emissions are slightly higher. Including the effect of aggregates caused a 40 % decrease in the simulated annual emissions relative to the simulations without accounting for the effects of aggregates. The new parameterization significantly improved the model performance regarding simulation of observed daily N 2O fluxes;r2 and RMSE improved from 0.11 and 198 g N 2O-N ha−1 d−1 to 0.41 and 40 g N2O-N ha−1 d−1, respectively. Our analyses of the model results show that aggregates have a larger impact on the reduction than on the production of N 2O. Reduction of N2O is more sensitive to changes in the drivers Correspondence to: P. C. Stolk (petra.stolk@wur.nl) than production of N2O and is in that sense the key to understanding N2O emissions from denitrification. The effects of changing environmental conditions on reduction of N 2O relative to N2O production strongly depend on the NO 3 content of the soil. More anaerobic conditions have hardly any effect on the ratio of production to reduction if NO 3 is abundant, but will decrease this ratio if NO 3 is limiting. In the first case the emissions will increase, whereas in the second case the emissions will decrease. This study suggests that the current knowledge of the hydrological, biogeochemical and physical processes may be sufficient to understand the observed N 2O fluxes from a fertilized clayey peatland. Further research is needed to test how aggregates affect the N 2O fluxes from other soils or soils with different fertilization regimes.


Introduction
Nitrous oxide (N 2 O) is a strong greenhouse gas (e.g.Denman et al., 2007) and negatively affects atmospheric ozone (Ravishankara et al., 2009).Agricultural soils are the largest anthropogenic source for N 2 O, at a global scale (Denman et al., 2007), in Europe (Schulze et al., 2009), and in The Netherlands (Van der Maas et al., 2008).Agricultural peatlands are major sources of N 2 O (Velthof and Oenema, 1995).
Observations of N 2 O emissions typically show a high temporal variability with long periods of low background emissions and a few high peak emissions lasting a couple of days to several weeks.Despite their short duration, these peaks contribute a major part to the total annual N 2 O emission (e.g.Scheer et al., 2008;Yamulki et al., 1995).When modelling N 2 O emissions, accurate simulation of peak emissions is therefore very important (e.g.Stehfest and Muller, 2004;Lamers et al., 2007).This requires models with time steps P. C. Stolk et al.: Modelling the effect of aggregates on N 2 O emission as small as one day or less.Various process-based field-scale simulation models for daily N 2 O fluxes are available (Chen et al., 2008).Whereas cumulative emissions have been simulated fairly well (e.g.Li, 2000;Saggar et al., 2007;Jarecki et al., 2008) the simulation of daily emissions is still poor (Groffman et al., 2009).
In a recent study with the SWAP-ANIMO model combination simulated N 2 O peak emissions from peatland occurred too early and were too high, compared to the observed N 2 O emissions (Stolk et al., 2011).It was hypothesized that this misfit was due to a violation of the assumption of direct equilibrium between the N 2 O concentration in the gaseous phase and in the aqueous phase.In SWAP-ANIMO first the total soil concentration is determined, after which the concentrations in the gaseous and aqueous phase are calculated, assuming equilibrium.However, if the assumption of direct equilibrium is not valid, the total soil concentration is not representative for the aqueous and gaseous concentration individually.In particular the gaseous concentrations need to be simulated realistically, since in unsaturated soils diffusion in the gas phase is the main mechanism for N 2 O transport to the surface.In other words, in cases of disequilibrium the total soil concentration cannot be used to derive N 2 O emission if the commonly applied equilibrium assumption is invoked.Support for this hypothesis is found in studies that prove a discrepancy between observed emissions and emissions calculated from soil concentrations using Fick's law for N 2 O (Neftel et al., 2007) and CO 2 (Koehler et al., 2010).
The assumption of direct equilibrium is commonly applied in biogeochemical models, although it is not always applicable in structured soils like clay or peat.Structured soils can contain aggregates with smaller pores and higher tortuosity, as well as disconnected or dead-end pores (Hoag and Price, 1997).These smaller or disconnected pores remain water filled longer than the larger pores in between the aggregates.The absence of air filled pores hampers the diffusion of oxygen (O 2 ) and results in formation of anaerobic microsites in the soil (Currie, 1961), favourable for denitrification (Groffman et al., 2009).Diffusion of N 2 O, an intermediate product of denitrification, out of these microsites is hampered for the same reason, increasing the residence time of N 2 O in the aqueous phase.This prolonged residence time facilitates further reduction of N 2 O to nitrogen gas (N 2 ).Thus, aggregates can both delay and reduce N 2 O emissions.
Various concepts have been proposed to describe soil structure and its effects on aeration and denitrification in simulation models, like aggregate models (Arah and Smith, 1989), pore models (e.g.Groenendijk et al., 2005), fractals (Rappoldt and Crawford, 1999), the multi-domain approach (Köhne et al., 2009), or the mobile-immobile approach (Van Genuchten and Wierenga, 1976).All of these model concepts have been used to simulate the effect of aggregates on nitrate (NO 3 ) and O 2 , but to the best of our knowledge none of them has been applied to simulate the effect of aggregates on N 2 O. Single-aggregate models, where the model-domain is limited to one aggregate only, have been used to study the effects of aggregates on the N 2 O concentration in the interpore space (Leffelaar, 1988;McConnaughey and Bouldin, 1985), but this work has not been extended to models that describe the soil as a collection of aggregates.
In the current N 2 O simulation models several concepts are applied that indirectly account for the effects of aggregates on N 2 O reduction and transport.For example, in the DNDC model-family (Li et al., 1992(Li et al., , 2000;;Li, 2000) a reduction factor decreases diffusion of N 2 O produced by denitrification based on the anaerobic fraction and the clay content.In PaSim (Schmid et al., 2001) N 2 O from denitrification is first assigned to a separate pool and a resistance is added for N 2 O from denitrification to enter the total soil N 2 O pool.The electron affinity for N 2 O relative to NO 3 was increased empirically to get more reduction of N 2 O to N 2 .The disadvantage of these concepts is that they do not describe the soil structure itself.Therefore each new situation may require a new calibration of the parameters involved in these concepts.Furthermore, effects of changes in soil structure are difficult to assess.
In this paper we present a parameterization to simulate the effects of aggregates on daily N 2 O emissions, following the mobile-immobile concept.In this approach the soil is separated into an immobile and a mobile zone, representing the pore space within and in between the aggregates, respectively.Vertical transport in the immobile zone is completely blocked, and occurs exclusively in the mobile zone.However, there is mass transfer between the zones.The transport between the zones is based on the shape and the size of the aggregates (Van Genuchten and Wierenga, 1976).The main advantage of the mobile-immobile approach is that the parameters are based on aggregate shape and aggregate size and can be derived from observations.The new parameterization is implemented in SWAP-ANIMO, a combination of a field-scale hydrological and biogeochemical model.
Simulations with the revised model are performed to test the hypothesis that aggregates cause delayed as well as reduced peak emissions and that the mobile-immobile model retrieves better simulation results than the original concept (Stolk et al., 2011) for a structured soil.To that end, the model simulations with the mobile-immobile model to account for the effect of aggregates are compared to simulations obtained with the original concept (Stolk et al., 2011) and with observations.In addition, a sensitivity study is carried out to determine how the aggregate size and shape affect production, reduction and emission of N 2 O and which model parameters affect these processes most.Finally we make a first assessment of the applicability and limitations of the new model concept for simulating N 2 O emissions.

General description of SWAP-ANIMO
SWAP (Soil-Water-Plant-Atmosphere) (Van Dam, 2000;Van Dam et al., 2008;Kroes et al., 2008) is a multi-layered simulation model for the transport of water, solutes and heat in unsaturated and saturated soils.The model is designed to simulate flow and transport processes at field scale level.Top, bottom, and lateral boundary conditions in SWAP allow computation of plant transpiration and interception, soil evaporation, runoff, irrigation, lateral drainage to and infiltration from drains and surface water, and seepage to or infiltration from deeper aquifers.
ANIMO (Agricultural NItrogen MOdel) (Rijtema et al., 1999;Groenendijk et al., 2005;Renaud et al., 2005;Hendriks et al., 2011) is a dynamic process-based simulation model for nutrients (N and P) and organic matter dynamics in the soil.In ANIMO vegetation interacts with soil moisture and nutrients.Aeration status is calculated from the average pore size, based on soil moisture content.The greenhouse gas module in ANIMO provides the simulation of production and consumption of N 2 O, CH 4 and CO 2 in the soil, as well as the vertical transport of these gases and their emission to the atmosphere.Transport of greenhouse gases involves molecular diffusion as well as advective transport by air and water flow.The minimum time step for ANIMO is one day; this time step is commonly used in greenhouse gas studies.Generally, a daily time step is sufficiently short to simulate N 2 O emissions with reasonable accuracy, because emission peaks generally last a few days to several weeks.Special events, like daytime-only thawing events, can cause peak emissions shorter than a day.In general, ANIMO is not suited to simulate those short-lived peak emissions.However, in ANIMO a special approach is implemented to account for short-lived peak emissions during and after precipitation events.In this approach for each precipitation event the fraction of the day the topsoil is saturated is determined.This fraction and the spatial anaerobic fraction are used to determine the overall aeration status of the top layer for that time step.This is explained in more detail in Appendix A.
ANIMO is coupled off-line to SWAP; the output from the latter model is used to prescribe water flow, soil moisture and soil temperature on a daily basis in ANIMO.Although the SWAP-ANIMO model combination has over 100 input parameters, for most soils in The Netherlands default values are available (Wolf et al., 2003;Wösten et al., 1994), which makes it easily applicable to soils in The Netherlands, and to comparable soils elsewhere, without extensive measurement campaigns or tuning of input parameters.

Original concept for N 2 O
The original ANIMO concept assumes instantaneous equilibrium between the air phase and water phase concentration of N 2 O throughout the soil.The ANIMO formulation of the conservation and transport (CT) equation for N 2 O reads (Hendriks et al., 2011): where and Here c s (kg N m −3 s ) and c w (kg N m −3 w ) are the N 2 O-N concentration in soil and soil water, respectively, t (d) is time, z (m) is depth, D eff (m 2 d −1 ) is the effective diffusion coefficient for simultaneous diffusion in soil air and water, q a (m 3 a m −2 d −1 ) and q w (m 3 w m −2 d −1 ) are the vertical air and water flux in the soil, respectively, α R (m 3 w m −3 a ) is the reciprocal of Bunsen's solubility coefficient, Q dr (m 3 w m −3 s d −1 ) is the discharge per layer to drainage systems and surface water, R pr,nit and R pr,den (kg N m −3 s d −1 ) are the N 2 O production rate by nitrification and denitrification, respectively, R rd (kg N m −3 s d −1 ) is the N 2 O reduction rate, θ a (m 3 a m −3 s ) and θ w (m 3 w m −3 s ) are the volumetric air and water content, respectively, p 1 and p 2 (-) are coefficients for vertical gaseous diffusion in soil, D 0,a and D 0,w (m 2 d −1 ) are the N 2 O diffusion coefficient in free air and free water, respectively, and τ w (-) is a tortuosity factor for vertical aqueous diffusion in soil.Emission to the atmosphere is composed of N 2 O diffusion over the atmosphere-soil boundary from soil water and soil air and advection of N 2 O by air flow.A more detailed description of N 2 O diffusion and emission, denitrification, nitrification and aeration in the original ANIMO is provided in Appendix A.

Mobile-immobile model concept
In the present study the mobile-immobile model concept has been implemented in SWAP-ANIMO for N 2 O. Figure 1 shows a schematic illustration of the mobile-immobile concept for N 2 O.The soil is divided into a mobile and an immobile soil fraction, F MO and F IM (-), respectively.All vertical transport is assumed to take place in the mobile zone.Nitrification is assumed to occur exclusively in the mobile soil fraction, whereas denitrification is assumed to occur exclusively in the immobile soil fraction.N 2 O may be transported from the immobile to the mobile zone, or vice versa, based on the concentration difference.The mobile-immobile approach is implemented for N 2 O only.The effects of soil structure on O 2 are already accounted for in the original ANIMO by a pore-model approach (Groenendijk et al., 2005).Furthermore, we assume that the effect of aggregates on NO 3 is negligible and that the NO 3 concentration is constant throughout the soil, because the diffusion coefficient for NO 3 in water is an order of magnitude higher than for N 2 O (Berner, 1971;Heincke and Kaupenjohann, 1999).We will come back to this assumption in the discussion.
The two main vertical transport processes for N 2 O, gaseous diffusive transport in soil air and convective flow with soil water, are leading in the determination of the size of the mobile and the immobile zone.Assuming that all air filled pores are interconnected, all air filled pore space is assigned to the mobile zone and the immobile zone is always water saturated.If structured soils are completely water saturated, part of the water is stagnant.This soil fraction with stagnant water represents the maximum immobile soil fraction, F IM,max (-).When the water filled pore space is larger than F IM,max , the immobile soil fraction, F IM , is equal to F IM,max : When water filled pore space is less than F IM,max , the immobile soil fraction is set to 95 % of the water filled soil fraction.This percentage denotes that always part of the water remains in the mobile soil fraction as water films around the particles (Quinton et al., 2009).Occurrence of anaerobic conditions required for denitrification is very unlikely when water filled pore space is less than the maximum immobile fraction.Consequently, this percentage does not significantly affect N 2 O emissions.
The total soil concentration c s is given by: Here c s,MO (kg N MO m −3 s ) and c s,IM (kg N IM m −3 s ) denote the concentrations of mobile and immobile N 2 O-N in the total soil, respectively.Whereas there can be disequilibrium between the N 2 O concentrations in the mobile and the immobile zone, within the mobile soil fraction we assume that the concentration in the soil air is in direct equilibrium with the local concentration in the soil water, so that we can write: From Eq. ( 6) we can derive the CT-equation for the N 2 O in the mobile zone: where Here D eff,MO (m 2 d −1 ) is the effective diffusion coefficient for the mobile zone and R tr (kg N m −3 s d −1 ) is the exchange rate of N 2 O between the mobile and the immobile zone.This term can be a source as well as a sink for the mobile zone.The 2nd and 3rd terms at the right hand side of Eq. ( 8) for convective transport and the 4th term for drainage have not been changed compared to the original CT-equation in Eq. (1).In these terms soil heterogeneity has already been taken into account in the soil specific hydraulic characteristics and therefore it is reasonable to assign all convective and drainage flow to the mobile zone.Also the nitrification term has not been changed compared to Eq. ( 1), because nitrification, which requires oxygen rich conditions, is likely to occur mainly in the mobile zone.
The CT-equation for N 2 O in the immobile zone includes only production and reduction of N 2 O by denitrification and exchange with the mobile zone: The exchange rate R tr in Eq. ( 10) is the same as in Eq. (8).

N 2 O production and reduction
Denitrification is the sequential reduction of nitrate (NO 3 ) to nitrogen gas (N 2 ) with N 2 O as an intermediate.It requires anaerobic conditions, so it encompasses both production and reduction of N 2 O.In ANIMO the ratio between production and reduction is determined by E af (-), the electron affinity of N 2 O relative to the electron affinity of NO 3 , as well as by the concentrations of N 2 O and NO 3 in water, c w,N 2 O and c w,NO 3 (kg N m −3 w ), respectively (Hendriks et al., 2011).A more detailed description of denitrification in general in ANIMO is provided in Appendix A.
In the mobile-immobile model the production and reduction rates are calculated as: Here R pr,el (kmol e − m −3 s d −1 ) is the electron (e − ) production rate, and 4/14 and 1/14 (kmol e − kg −1 N) express the electron equivalent ratio per molar mass of NO 3 -N and N 2 O-N, respectively.Functions f pH and f T (-) are response functions for pH and temperature (Hendriks et al., 2011).The NO 3 concentration is assumed to be equal in the mobile and immobile zones and all reduction of NO 3 is assumed to be located in the immobile zone.With respect to the original ANIMO we left out the response function for the effect of aeration on N 2 O production and reduction.This response function reduces the ratio of N 2 :N 2 O production if oxygen is present (Hendriks et al., 2011).In the current mobileimmobile model we assume that most oxygen that is still present in the soil is located in the mobile zone.The reduction function, which is based on the average oxygen content of the total soil, cannot be used for the immobile zone only.With regard to denitrification two situations can be distinguished: organic matter-limited or NO 3 -limited.If the organic matter content is limiting denitrification, the electron production rate, R pr,el , is known from organic matter decomposition and nitrification.If NO 3 is limiting, the production of N 2 O-N is equal to the reduction of NO 3 -N, which is calculated using a first order rate constant, k rd,NO 3 (d −1 ) (Hendriks et al., 2011).

Numerical elaboration
The CT-equations for the mobile and immobile zone are solved separately in an iterative approach.In the first iteration the transfer rate is calculated from estimates of the new concentrations in the mobile and immobile zone.Then alternately the CT-equations for the mobile and immobile zone are solved.Each iteration step the transfer rate is adjusted to the newly determined concentrations.The iteration is continued until the concentrations do not change significantly anymore.
The CT-equation for the mobile zone is solved for daily average concentrations.Accordingly, the CT-equation for the immobile zone is also solved for daily average concentrations.
The transfer rate of N 2 O from the immobile to the mobile zone is assumed to be constant over a day and is approximated as where k tr (d −1 ) is a mass transfer coefficient and θ w, * (m 3 w, * m −3 s ) is the moisture fraction of the zone with the highest concentration: θ w,IM if c w,IM > c w,MO and θ w,MO if c w,MO > c w,IM .This equation is an approximation of diffusive transport in soil water and has often been used in mobile-immobile zone studies, like Gerke and Van Genuchten (1993).For the description of the mass transfer coefficient k tr we follow Gerke and Van Genuchten (1993): Here β (-) is a shape factor depending on the geometry of the aggregates, a (m) represents the distance from the center of the aggregate to the interface with the mobile zone.
Assuming linear change with time during a timestep, t (d), the average concentration in the immobile zone can be approximated from Eq. ( 10) as Here c 0,w,IM is the concentration in the immobile zone at the start of the time step.Combining Eqs. ( 13) and ( 15) gives for the mass transfer rate: Further elaboration of Eq. ( 16) gives: with the mass transfer help variable φ (m 3 w,IM m −3 s d −1 ): This can be implemented in the CT-equation for the mobile zone (Eq.8): where The CT-equation for the mobile zone (Eq.19) is solved numerically with a tridiagonal system of equations, with the TRIDAG routine (Press et al., 1989).The same routine is used to calculate the concentrations at the end of the time step.The concentrations in the immobile zone at the end of the time step are calculated with Eq. ( 15), using t instead of 0.5 t.
Mass balance calculations to check for errors in the numerical implementation show no abnormalities.

New input parameters in the mobile-immobile model
The mobile-immobile model has three new input parameters to describe the soil structure: the maximum immobile soil fraction, F im,max , shape factor β and the half width of the aggregate a.Here we describe how these parameters can be determined and ranges reported in literature.
The maximum immobile soil fraction, F im,max , can be considered equal to the maximum fraction of stagnant water.This fraction can be determined experimentally with nuclear magnetic resonance measurements (Culligan et al., 2001), curve fitting of measured soil water content (Price and Whittington, 2010), 2-D image analysis with resin impregnation (Hoag and Price, 1997), or 3-D X-ray tomographic image analysis (Quinton et al., 2009).In these studies the maximum immobile soil fraction for peat was determined, ranging from 0.22-0.84.In a methane simulation model a similar parameter is used to define the immobile pore fraction for gaseous diffusion in peat, with a default value of 0.5 for non-tropical peatland (Walter and Heimann, 2000).For the immobile soil fraction in mineral soils, values in literature range from 0.04-0.98(Jaynes et al., 1995).
The value for shape factor β (-) depends on the type of aggregate.Values are available for a range of aggregate types and dimensions (Van Genuchten, 1985;Van Genuchten and Dalton, 1986).Typical values range from 3 for plane-sheet type aggregates to 15 for spherical aggregates.Cubic aggregates have a shape factor comparable to spherical aggregates.Cylindrical or prismatic aggregates have shape factors around 11. Larger shape factors are calculated for cylindrical or prismatic aggregates with a limited height.Smaller shape factors are calculated for hollow cylinders.
Parameter a (m) denotes the radius of a spherical aggregate, the half-width of a plane-sheet type aggregate, or the thickness of the soil cylinder surrounding a cylindrical pore.Typical sizes for fine to coarse aggregates in the FAO soil classification (FAO, 2006) are different for the different aggregate types.Typical values for β, a and β/a 2 (m −2 ) are  presented in Table 1, where β/a 2 represents a combined aggregate transfer factor.The larger this factor, the more rapid the transfer between the mobile and immobile zones.

Observations
Daily N 2 O emission data were available for the site of Oukoop.This is a grassland site in the western part of The Netherlands used for intensive dairy farming.The annual N application is about 350 kg N ha −1 and is a combination of manure and fertilizer.The water level in the ditches is kept at 45 to 55 cm below the average surface level.From August until November 2006, day of year (DOY) 230 to 310, 30 min.N 2 O fluxes have been measured with a combination of a quantum cascade laser spectrometer and a sonic anemometer following the eddy covariance method (Kroon et al., 2007).Daily average fluxes were calculated for days with >12 measurements.During this period fertilizer was applied once, on DOY 257.The topsoil consists of peaty clay and clayey peat on a subsoil of eutrophic peat.The soil was described in a pit and on four locations with augering, following the guidelines for soil profile description from the FAO (FAO, 2006).The soil was classified following the FAO classification system (FAO, 1998).The top layer (0-23 cm) has a fine to medium moderate granular and fine subangular  2 shows the precipitation, soil moisture content and fertilizer application during the observation period.

Values for input parameters
Values for the input parameters were equal to the values used in the former simulation study for Oukoop with the original ANIMO (Stolk et al., 2011).In that former study values for most input parameters were derived from observations, literature, existing databases, or the model defaults.The only parameters that were calibrated were two parameters for the hydraulic characteristics, the humus and organic matter content of the start-up run in 1941, two parameters describing nitrification, and E af , the electron affinity N 2 O relative to the electron affinity of NO 3 .
Values for the new input parameters β and a of the mobileimmobile model were based on the observations of the soil structure in the second soil layer.Shape factor β is set to 11 for all layers, representing a prismatic structure, with radius a of 50 mm.This set of parameters is thought to be representative for a wet soil, i.e. the conditions when denitrification occurs.F im,max was set to 0.5, the average value of the range reported in literature.

Analyses
The simulated N 2 O emissions with the mobile-immobile model for Oukoop were compared with the simulated emissions from the original model for Oukoop without aggregates.Simulated emissions were compared with the observed emissions and model performance was assessed with the generally used coefficient of determination r 2 (-), the RMSE (g N 2 O-N ha −1 d −1 ) (Neter et al., 1996), and the RMSE normalized to the standard deviation, RMSE n (-) (Kiese et al., 2005).Additionally, the modeling efficiency r 2 eff (-) (Nash and Sutcliffe, 1970) and the coefficient of residual mass, CRM (-), (Moreels et al., 2003) were determined.
In case of skewed distributions, as is often observed for daily N 2 O emissions, r 2 eff is a stricter measure than r 2 .A positive value (0 to 1) for r 2 eff indicates that the predicted value is a better estimate of the observations than the mean.CRM is the difference between the cumulative simulated and observed flux, relative to the cumulative observed flux and is used to assess the model performance to simulate cumulative fluxes.CRM is not affected by a time lag between observed and simulated fluxes.
A sensitivity analysis was performed for the Oukoop site to determine the most influential parameters for N 2 O production by denitrification, N 2 O reduction, and N 2 O emission and the difference in sensitivity of these processes between the topsoil (0-5 cm) and the subsoil (5-300 cm).The parameters used in the sensitivity analysis and their ranges are reported in Table 3.For each parameter the minimum and maximum values were determined based on literature, specific for structured soils (clay and peat).The base values are the input values of the Oukoop model.A reference run was performed with all parameters at their base value.In an One-At-a-Time-analysis (Campolongo et al., 2001) each parameter was set alternately at its minimum and maximum value, keeping all other parameters constant.For each run the cumulative N 2 O emission was determined, as well as the cumulative N 2 O production and reduction for the topsoil and the subsoil.The sensitivity is presented as the cumulative production, reduction and emission, relative to the cumulative production, reduction and emission of the reference run, respectively.The parameters β and a are taken together as β/a 2 .The effect of the new model parameters on the emissions was evaluated.

Comparison of model results
Figure 3 shows the emissions simulated for the Oukoop-site with and without taking into account the effect of aggregates.This figure clearly shows the improvement of the model's capability to simulate the observed fluxes.The effect of aggregates on the emission is as expected: smaller peak emissions due to increased reduction, a delay in most peaks, and higher emissions in between the peaks due to a longer storage time in the immobile zone.Generally the onset of the peak emissions has not changed.
Table 4 shows the statistics indicating the performance of both model concepts regarding the simulation of the observations.The positive value of r 2 eff (0.30) for the model with aggregates indicates that those simulation results are reasonably well, in spite of the occurrence of emission peaks which often causes very low or negative values of that indicator (cf. the simulation without aggregates).The simulated cumulative emissions for the simulation period is only slightly too high (CRM = 0.03).The simulated annual emission for 2006 (24.5 kg N 2 O-N ha −1 ) is in good agreement with the annual emission that has been derived from observations (17.8 kg N 2 O-N ha −1 ± 58 %) (Kroon et al., 2010).This is a more than 40 % decrease compared to the annual emission simulated with the model without aggregates (42.0 kg N 2 O-N ha −1 ) and clearly a major improvement as well.
The net effect of the aggregates on the emissions can be separated into two components: (1) a decrease in the peak emissions due to enhanced reduction of N 2 O and (2) an increase of the "background emissions" in between the peaks due to a longer storage time.Both components contribute to the difference between the annual emission with and without aggregates (17.4 kg N 2 O-N ha −1 ) although with an opposing effect: the annual cumulative peak emissions decrease with 22.2 kg N 2 O-N ha −1 , whereas the annual cumulative background emissions increase with 4.7 kg N 2 O-N ha −1 .Total gaseous N production is larger in the aggregate model, whereas NO 3 reduction is less.The fraction of N 2 O/(N 2 O + N 2 ) is smaller in the aggregate model, which is illustrated in Fig. 4 for the topsoil.
In the period between DOY 250 and 270 the simulated flux is relatively low compared to the observed flux, both with and without aggregates.During the precipitation period from DOY 230 to 246 production mainly is simulated in the deeper soil layers.Most of the N 2 O produced there initially dissolves in the soil water, because in these layers the high (>0.99)water filled pore space, θ w /θ sat , hampers gaseous diffusion.The dissolved N 2 O is then reduced further to N 2 .In the aggregate model the emission is slightly higher, because in the mobile zone reduction cannot take place.

Sensitivity analysis
In Fig. 5a-c the simulated daily N 2 O emissions with effect of aggregates are plotted for various values for the new model parameters F IM,max , β and a, respectively.A smaller immobile soil fraction (Fig. 5a) in general gives smaller emissions.Changes in F IM,max do hardly affect N 2 O production by denitrification, but they do affect reduction.When the immobile soil fraction is smaller and the N 2 O production in the immobile soil fraction is the same, the N 2 O concentration in the immobile zone will be higher.This results in more reduction of N 2 O to N 2 .
Less mass transfer between the mobile and immobile zone causes smaller peak emissions.In this case, less mass transfer can result from either a smaller shape factor (Fig. 5b) or a larger radius (Fig. 5c).The start of the peak emissions does not change, but timing of the highest peaks might occur one day later when there is less mass transfer.The emissions in between the peaks are higher when mass transfer is less.Decrease of the peak emissions is caused by increased N 2 O reduction; N 2 O production is hardly affected.
Figure 6 clearly shows that in the simulations aggregate size has a large impact on N 2 O reduction.In soils with large aggregates there is more reduction than in soils with small aggregates.This in turn affects the emission; more reduction leads to less N 2 O emission from soils with large aggregates, compared to soils with small aggregates (Fig. 5c).In soils with small aggregates, or in the absence of aggre- gates, peak emissions are therefore much higher.The difference in emissions from soils with a different shape factor can also be traced back to a difference in the reduction of N 2 O (not shown).For soils with aggregates with a radius >2 cm, the changes in the simulated emissions are considerable and therefore it is necessary to take aggregates into account in the simulation of N 2 O emissions.This is in line with Arah (1990), who stated that simulation of denitrification in soils with aggregates >1 cm requires a model that takes spatial heterogeneity into account.
A few experiments have considered the effect of aggregate size on denitrification (Miller et al., 2009;Drury et al., 2004;Seech and Beauchamp, 1988).The results were contradicting, with more denitrification, less denitrification or no significant difference in denitrification with increasing aggregate size.However, in all experiments N 2 O reduction to N 2 was blocked.So, there is no experimental evidence yet to validate our conclusions on the effect of aggregates size on N 2 O reduction and emission.
The results of the general sensitivity analysis are presented in Fig. 7.The electron affinity E af has the largest impact on the emission, followed by the aggregate transfer factor β/a 2 .The other parameters have less effect on the emissions.The electron affinity and the mass transfer factor, and to a lesser extent F IM,max , have hardly any impact on the production of N 2 O, but a large impact on the reduction of N 2 O, mainly in the topsoil.The coefficient for gaseous diffusion in the top layer p 2 (1) affects both production and reduction of N 2 O in the topsoil and the subsoil.Interestingly, the effect in the topsoil is opposite to the effect in the subsoil.The rate constants for nitrification, k ox,NH 4 and denitrification, k rd,NO 3 , also affect both production and reduction of N 2 O.The rate constant for denitrification has most effect on reduction on the topsoil, whereas the rate constant for nitrification has the least effect on production of N 2 O in the topsoil.The hydraulic conductivity of the top layer, K sat (1), has a relatively large effect on reduction and production in the topsoil.
It is interesting to notice that reduction of N 2 O is more sensitive to changes in the parameters than N 2 O production, and has more influence on the emissions.reduction of N 2 O to N 2 is the key to understanding N 2 O emissions from denitrification.
It is interesting to note that reduction in the topsoil is more sensitive to parameter changes than reduction in the subsoil.As pointed out before, in the subsoil typically the water filled pore space >0.99 during denitrification.Here, anaerobicity is mainly caused by limited diffusion of oxygen in the nearly saturated soil, which also hampers diffusion of N 2 O. Therefore, both with and without aggregates, most N 2 O that is produced in deeper soil layers is reduced further to N 2 .In the topsoil typically the water filled pore space <0.99 during denitrification.Here, anaerobicity is mainly caused by temporal saturation of the soil, during and shortly after strong precipitation events.After the excess precipitation has infiltrated, gaseous diffusion of N 2 O restarts immediately.Because of the "competition" for N 2 O between emission and reduction, changes in reduction related to aggregates have more effect in the topsoil than in the subsoil.
Another interesting aspect in Fig. 7 is that both the minimum and the maximum value of K sat (1) result in lower emissions than the base value.K sat is the saturated conductivity and in the top layer it determines the infiltration rate at the soil surface.A low K sat hampers infiltration, which enhances saturated conditions.As described before, this results in more anaerobicity in the top layer.For high values of K sat the opposite occurs.A more detailed description of this process is provided in Appendix A.
The effect of the changing aeration status of the topsoil on N 2 O production and reduction is further explored in Fig. 8.This figure shows the results of simulations with a low and high saturated conductivity in the top layer (Fig. 8a and b, respectively), representing more and less anaerobicity in the top layer.It can be seen that the more anaerobic conditions in the topsoil (Fig. 8a) can cause both an increase and a decrease in the peak emissions compared to the less anaerobic conditions (Fig. 8b).This is related to the two situations that can be distinguished for denitrification: organic matter limited or NO 3 limited.Around DOY 283 and 291 more anaerobicity leads to an increase in peak emissions.In these cases NO 3 is not limiting and both NO 3 reduction and N 2 O reduction proceed at a rate higher than in the situation with less anaerobicity (Fig. 8b).The ratio between NO 3 and N 2 O reduction is hardly affected.Around DOY 275, 278, and 296 more anaerobicity leads to a decrease in peak emissions.In these cases the available NO 3 is insufficient to meet the potential organic matter decomposition.Consequently, the ratio between NO 3 and N 2 O reduction is affected, in favour of N 2 O reduction.So, the effect of anaerobicity on the ratio between NO 3 and N 2 O reduction is opposite for situations with NO 3 excess or NO 3 shortage.Recently, this contrasting effect of NO 3 content on the ratio of N 2 O production and reduction was reported in a field study for a wet and a dry section of a natural peatland (Roobroeck et al., 2010).Figure 8 also reveals that changes in K sat (1) can affect the onset of peak emissions, for example around DOY 275.In the case with the higher conductivity (Fig. 8b) this improves the model performance considerably, compared to the default simulation with aggregates (r 2 eff = 0.50 versus 0.30 for the default case).Interestingly, the daily soil moisture content in the top soil has hardly changed.This shows the importance of short periods of saturation of the top soil during and after precipitation events for N 2 O emissions.

Applicability and limitations
In this study we have assumed that the effect of aggregates on NO 3 is negligible and that the NO 3 concentration is constant throughout the soil.The good model performance does not cast a doubt on this assumption.However, other model studies have shown that aggregates do affect the NO 3 concentrations; NO 3 concentrations within the aggregates are lower than in between them (Leffelaar, 1988;Arah, 1990).As we discussed in the foregoing, our analyses show that NO 3 content can have a large effect on the ratio between N 2 O production and N 2 O reduction: in case of NO 3 excess this ratio is relatively stable, whereas in case of a shortage of NO 3 this ratio can change.Especially in unfertilized soils with little NO 3 , a difference between the NO 3 concentration within and in between the aggregates might cause a major shift in the ratio between N 2 O production and N 2 O reduction and thus in the N 2 O emission.In these cases it is necessary to include NO 3 in the mobile-immobile model.The results of the present study suggest that in fertilized soils the effect of aggregates on the NO 3 concentration hardly affects the N 2 O emissions and that the current parameterization is adequate.

Conclusions
The results of this study strongly confirm that aggregates affect simulated N 2 O emissions.More in general, this study emphasizes the importance of physical transport processes in N 2 O simulation models.The main effect of the aggregates is to increase the reduction of N 2 O to N 2 .This reduces the strength of the peak emissions.Longer storage of N 2 O also leads to slightly higher background emissions.Reduc-tion of N 2 O is more variable than production of N 2 O and is in that sense the key to understanding N 2 O emissions from denitrification.According to our simulations, the NO 3 content has a major impact on the ratio between N 2 O production and N 2 O reduction and how this ratio is affected by changing environmental conditions.
Implementation of the effect of aggregates in the SWAP-ANIMO model combination highly improved the model performance to simulate the observed N 2 O fluxes from a fertilized peat soil.With our implementation of these processes in the SWAP-ANIMO model combination we were able to simulate the observed fluxes reasonably well without extensive model calibration.As such this model combination is a promising tool for inventories, scenario studies or more integral scientific studies.However, model evaluation so far has been limited to data series of one site only.As a next step the model should be validated for other sites with different soil types, water management, fertilizer application, or climate.

Fig. 1 .
Fig. 1.Conceptual representation of the mobile-immobile model concept for N 2 O.
MO c w,MO and c s,IM = θ w,IM c w,IM (6) where θ w,MO = θ w − θ w,IM and θ w,IM = F IM θ sat (7) Here θ w,IM (m 3 w,IM m −3 s ) and θ w,MO (m 3 w,MO m −3 s ) represent the immobile and mobile moisture content relative to the total soil, respectively, and c w,IM (kg N m 3 w,IM ) and c w,MO (kg N m 3 w,MO ) represent the concentrations of N 2 O-N in water in the mobile and immobile zone, respectively.

Fig. 4 .
Fig. 4. Simulated fraction N 2 O of total N 2 + N 2 O production against water filled pore space (wfps) in the topsoil.

Fig. 5 .
Fig. 5. N 2 O emission for a range of values of (a) the maximum immobile fraction F im,max (= fim); (b) the shape factor β (= bet); and (c) the radius a (= alf) of a cylindrical aggregate.All other parameters are at their base value.

Fig. 6 .
Fig. 6.N 2 O reduction in the topsoil (0-5 cm) for a range of values for the radius a (= alf).

Fig. 7 .
Fig. 7.Results of the sensitivity analysis: range of cumulative N 2 O emission, production by denitrification, and reduction, relative to the cumulative N 2 O emission, production, or reduction of the reference run.The ranges for production and reduction are given separately for the topsoil (0-5 cm) and the subsoil (5-300 cm).The open and filled circles represent the minimum and maximum value of the parameter, respectively.

Table 1 .
Typical values for the input parameters β and a. a denotes the radius of a spherical or cylindrical aggregate, or the half-width of a plane-sheet type aggregate. *

Table 3 .
Ranges and base values for the main N 2 O parameters used in the sensitivity analysis.

Table 4 .
Model performance regarding the simulation of the observed daily N 2 O emissions.

Table A1 .
Description of the symbols used in Appendix A. Q 10,p/ox,NO 2 ratio Q 10 factor NO 2 production to Q 10 factor NO 2 oxidation 2 -Q 10,rd,N 2 O/NO 3 ratio Q 10 factor N 2 O reduction to Q 10 factor NO 3 reduction −3 d −1 R pr,el,an electron production rate in anaerobic soil fraction kmol e − m −3 d −1 R pr,el,an,0 potential electron production rate in anaerobic soil fraction kmol e − m −3d −1 R pr,N 2 O,den production rate N 2 O during denitrification kg N m −3 d −1 R pr,N 2 O,nit production rate N 2 O during nitrification kg N m −3 d −1 R rd,N 2 O reduction rate N 2 O kg N m −3 d −1 R rd,O 2oxygen reduction rate from decomposition and nitrification kg O m −3 d −1