To what extent can soil moisture and soil Cu contamination stresses affect nitrous species emissions? Estimation through calibration of a nitrification–denitrification model

Continental biogeochemical models are commonly used to predict the effect of land use, exogenous organic matter input or climate change on soil greenhouse gas emission. However, they cannot be used for this purpose to investigate the effect of soil contamination, while contamination affects several soil processes and concerns a large fraction of land surface. For that, in this study we implemented a commonly used model estimating soil nitrogen (N) emissions, the DeNitrification DeCompostion (DNDC) model, with a function taking into account soil copper (Cu) contamination in nitrate production control. Then, we aimed at using this model to predict N₂O-N, NO₂-N, NO-N and NH₄-N emissions in the presence of contamination and in the context of changes in precipitations. Initial incubations of soils were performed at different soil moisture levels in order to mimic expected rainfall patterns during the next decades and in particular drought and excess of water. Then, a bioassay was used in the absence or presence of Cu to assess the effect of the single (moisture) or double stress (moisture and Cu) on soil nitrate production. Data of nitrate production obtained through a gradient of Cu under each initial moisture incubation were used to parameterise the DNDC model and to estimate soil N emission considering the various effects of Cu. Whatever the initial moisture incubation, experimental results showed a NO₃-N decreasing production when Cu was added but depending on soil moisture. The DNDC-Cu version we proposed was able to reproduce these observed Cu effects on soil nitrate concentration with r² > 0.99 and RMSE < 10 % for all treatments in the DNDC-Cu calibration range (> 40 % of the water holding capacity) but showed poor performances for the dry treatments. We modelled a Cu effect inducing an increase in NH₄-N soil concentration and emissions due to a reduced nitrification activity and therefore a decrease in NO₃-N, N₂O-N and NO_x-N concentrations and emissions. The effect of added Cu predicted by the model was larger on N₂-N and N₂O-N emissions than on the other N species and larger for the soils incubated under constant than variable moisture. Our work shows that soil contamination can be considered in continental biogeochemical models to better predict soil greenhouse gas emissions.

1 Introduction

The increase in atmospheric greenhouse gases (GHG) like CO₂, CH₄ or N₂O is expected to induce a global climate change with higher mean temperature or changes in rainfall patterns with projections of increased precipitations or droughts depending on regions (Knutti and Sedláček, 2012). These modifications in rainfall patterns may impact soil moisture, which is one of the main drivers of soil microbial activity (Moyano et al., 2013; Schimel, 2018; Stark and Firestone, 1995). Microbial communities ensure key activities supporting numerous ecosystem functions, such as those involved in nitrogen (N) cycle influencing N₂O emissions (Butterbach-Bahl et al., 2013; Galloway et al., 2008) and are at the origin of more than 80 % of N₂O fluxes (IPCC, 2019). In particular, nitrification–denitrification processes are largely controlled by the local (an-)oxic treatments and therefore by soil moisture (Borken and Matzner, 2009; Fierer et al., 2003; Guo et al., 2014; Schimel, 2018), denitrification being the main source of soil N₂O emission for moist soils whereas for dry soils N₂O emissions are mainly due to nitrification (Bateman and Baggs, 2005). N soil flux dynamics are thus particularly difficult to predict at a large scale because of this strong dependency on local soil O₂ availability (Khalil et al., 2004). Despite this, some continental biogeochemical models have shown improved predictions when the N cycle was explicitly represented (Butterbach-Bahl et al., 2009; Kesik et al., 2005; Vuichard et al., 2019).

In addition to climate change, human activities introduce significant quantities of contaminants into the environment, such as trace elements (TEs), which are persistent and can be toxic for soil biota (Bech et al., 1997; Giller et al., 2009). Indeed, the contamination of soils by TE has become a major concern at the global scale (De Vleeschouwer et al. 2007; Khan et al. 2008) coming from atmospheric sources (Steinnes et al., 1997) or through the use of pesticides (Nicholson et al., 2003). In particular, TE contaminations are known to largely affect soil microorganisms (Bååth, 1989; Giller et al., 2009) and their activities, such as nitrification–denitrification processes (Broos et al., 2007; Mertens et al., 2010). Therefore, the combined effect of climate change and of soil contamination may largely impact the emissions of NO_x and N₂O from soils (Holtan-Hartwig et al., 2002; Vásquez-Murrieta et al. 2006). However, the effect of the interactions between climate change and soil contamination on the GHG emissions is still poorly documented (Rillig et al., 2019; Zandalinas et al., 2021).

Despite recent progress, the Earth system models (ESMs) used to predict future climate change still do not take into account the soil contamination effect on GHG emissions (Anav et al., 2013) in spite of the fact that at a large spatial scale many soils are listed as contaminated (Rodríguez-Eugenio, 2018; Lado et al., 2008). Furthermore, soil biogeochemical models are often used to estimate loss or accumulation of N species (ammoniac NH₄ volatilisation, nitrate NO₃ leaching) (Giltrap et al., 2010) or their respective concentrations under scenarios of organic fertiliser amendments but do not take into account the contamination which often occurs simultaneously (Wuana and Okieimen, 2011). Thus, there is a growing need to provide continental models combining ecotoxicological/contamination and climate change concerns. Among the biogeochemical models, DeNitrification DeCompostion (DNDC, Li et al., 1992) is a relatively simple model handling both biogeochemistry of denitrification and microbial growth (Li et al., 2000), and on which the land surface model soil N component – a part of ESMs like ORCHIDEE – is built (Vuichard et al., 2019).

In order to improve model outputs, this study combines, in an innovative way, experimental and modelling approaches to evaluate the impact of soil moisture on the sensitivity of nitrification to copper (Cu) toxicity and consequently on GHG N emissions. Cu was chosen as a model of soil contamination due to both its relevance in agricultural soils and available data in the literature (Broos et al., 2007; Mertens et al., 2010; Sauvé et al., 1999). It is not straightforward to assess distinct effects between punctual or chronic contamination on microbial structure or soil functions (Brandt et al., 2010; Oorts et al., 2006; Smolders et al., 2009). Here, we designed experiments to assess the conjugated effects of trace metal contamination and soil moisture stress on soil N cycle. Soil initial incubations were run for 5 weeks by applying a given soil moisture from drought to water saturation. Then, a bioassay with a gradient of Cu added by spiking was performed to estimate ${\mathrm{NO}}_{\mathrm{3}}^{-}$ production. The experimental data were used to calibrate a new model, DNDC-Cu, able to predict NO_x and N₂O emissions with the implementation of new functions considering the effect of Cu concentration ([Cu]) on nitrification–denitrification processes. Our hypothesis is that the building of such a model allows a gain in the understanding of the effect of a soil [Cu] on NO_x and N₂O and NH₄ cycling in a climate change context. Hence, data are also used here to discuss knowledge gaps in such modelling approaches and to question the matter of soil contamination data in climate change scenarios.

Figure 1 Schematic representation of the experimental and modelling procedures. Left refers to the experimental part and centre to right to the modelling part. Soils were first incubated 5 weeks at different constant percentages of the water holding capacity (WHC) or at two variable moisture levels, Dry-Only (DO) and Dry-Rewet (DR). Then NO₂-N, NO₃-N and NH₄-N soil concentrations were measured after this initial incubation, and values were used to initialise DNDC, while a bioassay was also applied to soil aliquots. The 3 d bioassay included ${\mathrm{NH}}_{\mathrm{4}}^{+}$ in excess and copper (Cu) spikes at 0, 50, 100, 250, 500, 750, 1000 and 2000 mg Cu kg soil⁻¹ of soil. After 1 and 3 d of bioassay incubation, NO₃-N production was measured in the supernatant. NO₃-N productions against [Cu] gradients were used to define the functions of Eqs. (28) to (31) in Sect. 3.1 (see text). Soil respiration values were extracted from the curve C_i of Fig. 1 in Annabi et al. (2007).

2 Materials and methods

2.1 Soil sampling

The soil was sampled in January 2017 at the surface layer (0–20 cm) of a control plot at the QualiAgro experimental site (48^∘87^′ N, 1^∘97^′ E – https://www6.inrae.fr/valor-pro_eng/Experimental-devices/QualiAgro/QualiAgro-web-site, last access: 2 November 2021). The soil sample was immediately wet sieved at 5 mm and shortly stored at 4 ^∘C until microcosm build-up. Aliquots of this sieved soil were used to measure the initial water content in addition to the maximum water holding capacity (WHC) for the further microcosm experiments. This site is located at Feucherolles near Paris, France, and had been designed to evaluate urban compost fertility together with the monitoring of contaminant inputs (Cambier et al., 2019). Soil is a Luvisol with 15 % clay, 78 % silt and 7 % sand; a pH of 6.9; organic carbon (Corg) and total N contents at 10.5 ± 0.2 and 1.00 ± 0.03 g kg⁻¹ soil, respectively; and a catatonic exchange capacity (CEC) of 7.9 ± 0.8 ${\mathrm{cmol}}^{+}{\mathrm{kg}}^{-\mathrm{1}}$ soil. This soil is not contaminated with Cu, and geochemical [Cu] background measured by inductively coupled plasma atomic emission spectroscopy (ICP-AES) after HF-HClO₄ extraction was 12 mg Cu kg⁻¹ soil.

2.2 Experimental setup

In order to evaluate the impact of soil moisture on the sensitivity of nitrification to Cu toxicity, we carried out a two-step experiment. The first step consisted in initial incubations at five different WHCs over 5 weeks, and the second step was a 3 d bioassay with a spiked Cu gradient (Fig. 1).

For the 5-week initial incubation, five microcosms were built up with about 5 g of sampled soil. Three of them were set up with a constant moisture corresponding to 30 %, 60 % and 90 % of their WHC in order to span, respectively, limiting, optimal and saturating conditions for the microbial activities. These three samples will hereafter be called 30 %, 60 % and 90 %, respectively. Their water contents were verified by weighting every 2 d and water was added if necessary. The two other microcosms were incubated in order to simulate two kinds of drought and dry–rewet cycles. One, hereafter called “Drought” (or DO), started with 1 week at 60 % WHC, and then the soil was left for 3 weeks without added water to mimic a dry period until 10 % of the WHC before rewetting at the initial 60 % WHC. The other, hereafter called “Dry-Rewet” (or DR) encountered two cycles of a 1-week near-saturation period (90 % WHC) followed by a 1-week dry period (10 % of the WHC) and ending with a 1-week near-saturation period. Drying was performed by natural evaporation (gentle air-drying at the laboratory temperature, i.e. 20 ^∘C) and controlled by weighting.

At the end of the initial incubation period, we performed a nitrification bioassay using three replicates originating from soils and following an adaptation of the method proposed by Petersen et al. (2012). The bioassay consisted in nitrate production measurement over a short-term aerobic incubation in soil slurries (ratio of soil to solution was 1:10) with ammonium in excess and in the presence of gradients of Cu. Briefly, 3.5 g of fresh soil (approximately 3 g of soil equivalent dry weight) was mixed in 50 mL Falcon^® tubes with 29 mL of a 10 mM HEPES buffer solution (hydroxyethyl piperazineethanesulfonic acid, Sigma-Aldrich, France) to maintain a constant pH under Cu spiking and nitrification activity and containing the substrate (NH₄)₂SO₄ (3 mM) (Sigma-Aldrich, France). Soils were first spiked with a gradient of increasing Cu²⁺ in the presence of an excess of ${\mathrm{NH}}_{\mathrm{4}}^{+}$, and the resulting potential nitrification activity (PNA) was measured.

The microcosms incubated at constant moisture were kept at their moisture level (30 %, 60 % or 90 % of WHC), whereas those incubated at variable moisture levels were set at 60 % WHC. The ${\mathrm{NO}}_{\mathrm{3}}^{-}$production rates were measured in soil slurries over a short-term aerobic incubation, for each Cu added concentration. Briefly, 1 mL of Cu solution at different concentrations was added in soil slurries to reach added [Cu] of 50, 100, 250, 500, 750, 1000 and 2000 mg Cu kg⁻¹ of soil (final soil [Cu] of 62, 112, 262, 512, 762, 1012 and 2012 mg Cu kg⁻¹ of soil and control with 12 mg Cu kg⁻¹ of soil). The pH was adjusted to 7. Then, microcosms were incubated on a rotary shaker (150 rpm) under aerobic conditions at 25 ^∘C until 72 h. After 0, 24 and 72 h of incubation, 2 mL aliquots of 3 g were transferred in Eppendorf vials and centrifuged. The supernatants were collected and stored in microplates at −20 ^∘C until analyses of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ by colorimetric determinations, following the reduction of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ in ${\mathrm{NO}}_{\mathrm{2}}^{-}$ by vanadium(III) and then the detection of ${\mathrm{NO}}_{\mathrm{2}}^{-}$ by the acidic Griess reaction (Miranda et al., 2001). Finally, PNA ($\mathrm{\mu }\mathrm{g}{\mathrm{NO}}_{\mathrm{3}}\text{-}\mathrm{N}{\mathrm{g}}^{-\mathrm{1}}\mathrm{soil}{\mathrm{h}}^{-\mathrm{1}}$) was calculated on the basis of ${\mathrm{NO}}_{\mathrm{3}}^{-}\text{-}\mathrm{N}$ + ${\mathrm{NO}}_{\mathrm{2}}^{-}\text{-}\mathrm{N}$ concentrations measured at different time steps. In our bioassay, [${\mathrm{NO}}_{\mathrm{2}}^{-}$] was negligible and PNA was thus calculated following Eq. (1), by checking the linear production rate of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ between 2, 24 and 72 h:

$$\begin{array}{}\text{(1)}& \text{PNA}={\displaystyle \frac{[{\mathrm{NO}}_{\mathrm{3}}^{-}{]}_{\text{Tfinal}}-[{\mathrm{NO}}_{\mathrm{3}}^{-}{]}_{\text{Tinitial}}}{\text{Tfinal}-\text{Tinitial}}}\times V_{\mathrm{s}}÷W,\end{array}$$

where V_s is volume of solution, W is weight of fresh soil and T is time of incubation.

Cu in solution was measured by centrifugation of the soil–solution mixture of each bioassay, followed by a determination of Cu in solution by flame atomic absorption spectroscopy. Cu-in-solution values are provided in Table S1 in the Supplement.

2.3 Nitrification–denitrification model

Nitrification and denitrification processes are represented following the DNDC model proposed by Li et al. (1992, 2000). In this study, we used a simplified version of DNDC adapted by Zaehle and Friend (2010) initially calibrated for soil WHC > 40 %, which we intended here to test for 30 % of WHC. This simplified version needs less boundary data but keeps a mechanistic description of the main processes. Modelled N species are expressed in amount of N, i.e. NH₄-N, NO₃-N, NO_x-N and N₂O-N. To be able to represent both nitrification and denitrification processes occurring in aerobic and anaerobic sites, the soil is split into aerobic and anaerobic fractions based on an empirical relationship linking O₂ consumption to soil respiration. In aerobic microsites, nitrification takes place following Eq. (2):

$$\begin{array}{}\text{(2)}& \begin{array}{rl}\text{nitrification}=& f\left(\text{SWC}\right)\times f\left(\text{temp}\right)\times f\left(\text{pH}\right)\times k_{\text{Nit}}\\ & \times (\mathrm{1}-\text{anvf})\times {\mathrm{NH}}_{\mathrm{4}},\end{array}\end{array}$$

with NH₄-N being the stock of ammonium (in g N m⁻²); (1−anvf) the aerobic fraction of the soil described thereafter in Eq. (21); k_Nit the nitrification rate (d⁻¹); f(SWC), f(temp) and f(pH) three rate modifiers representing the effect of soil water content (m³ m⁻³); and temperature (K) and pH as scalars, respectively. They are described by the following Eqs. (3)–(5).

$$\begin{array}{}\text{(3)}& {\displaystyle }\begin{array}{rl}f\left(\text{SWC}\right)=& \mathrm{0.0243}+\mathrm{0.9975}\times \text{SWC}+\mathrm{5.6358}\times {\text{SWC}}^{\mathrm{2}}\\ & +\mathrm{17.651}\times {\text{SWC}}^{\mathrm{3}}+\mathrm{12.904}\times {\text{SWC}}^{\mathrm{4}}\end{array}\text{(4)}& {\displaystyle }\begin{array}{rl}f\left(\text{temp}\right)=& \mathrm{0.0233}+\mathrm{0.3094}\times \text{temp}+\mathrm{0.2234}\times {\text{temp}}^{\mathrm{2}}\\ & +\mathrm{0.1566}\times {\text{temp}}^{\mathrm{3}}+\mathrm{0.0272}\times {\text{temp}}^{\mathrm{4}}\end{array}\text{(5)}& {\displaystyle }f\left(\text{pH}\right)=\mathrm{1.2314}+\mathrm{0.7347}\times \text{pH}+\mathrm{0.0604}\times {\text{pH}}^{\mathrm{2}}\end{array}$$

The NH₄-N nitrified is transformed into N₂O-N, NO-N or NO₃-N due to microbial processes and chemonitrification following Eqs. (6)–(8).

$$\begin{array}{}\text{(6)}& {\displaystyle }{\mathrm{\text{Nitrification}}}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}=\text{ftv}\times \text{SWC}\times k_{{\text{Nitrif}}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}}\times \text{nitrification}\text{(7)}& {\displaystyle }\begin{array}{rl}{\text{Nitrification}}_{\mathrm{NO}}=& \text{ftv}\times \text{SWC}\times k_{{\text{Nitrif}}_{\mathrm{NO}}}\times \text{nitrification}\\ & +\mathrm{496}\mathrm{950}\times e^{-\mathrm{1.62}\times \text{pH}}\\ & \times e^{-\mathrm{31}\mathrm{494}/(\text{temp}\times R)}\\ & \times \text{nitrification}\end{array}\text{(8)}& {\displaystyle }\begin{array}{rl}{\text{Nitrification}}_{{\mathrm{NO}}_{\mathrm{3}}}=& \text{nitrification}-{\text{nitrification}}_{\mathrm{NO}}\\ & -{\text{nitrification}}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}\end{array}\end{array}$$

Here $k_{{\text{Nitrif}}_{\mathrm{NO}}}$ and $k_{{\text{Nitrif}}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}}$ are two fixed rates (d⁻¹), ftv a rate modifier controlled by temperature and given in Eq. (9), and R the ideal gas constant.

$$\begin{array}{}\text{(9)}& \text{ftv}={\mathrm{2.72}}^{\left(\mathrm{34.6}-\frac{\mathrm{9615}}{\text{temp}}\right)}\end{array}$$

Then, the NO₃-N produced during the nitrification process enters the denitrification module where it is reduced sequentially into NO_x-N, N₂O-N or N₂-N following Eqs. (10) to (12).

$$\begin{array}{}\text{(10)}& {\displaystyle }& {\displaystyle }{\mathrm{\text{Denitrification}}}_{{\mathrm{NO}}_{\mathit{\text{x}}}}=\text{anvf}\times \left({\displaystyle \frac{{\mathit{\mu }}_{{\mathrm{NO}}_{\mathrm{3}}}}{\mathrm{0.401}}}+\mathrm{0.09}\times {\displaystyle \frac{{\mathrm{NO}}_{\mathrm{3}}}{{\mathrm{N}}_{\mathrm{tot}}}}\right)\times B\text{(11)}& {\displaystyle }& {\displaystyle }{\mathrm{\text{Denitrification}}}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}=\text{anvf}\times \left({\displaystyle \frac{{\mathit{\mu }}_{{\mathrm{NO}}_{\mathit{\text{x}}}}}{\mathrm{0.428}}}+\mathrm{0.035}\times {\displaystyle \frac{{\mathrm{NO}}_{\mathit{\text{x}}}}{{\mathrm{N}}_{\mathrm{tot}}}}\right)\times B\text{(12)}& {\displaystyle }& {\displaystyle }{\mathrm{\text{Denitrification}}}_{{\mathrm{N}}_{\mathrm{2}}}=\text{anvf}\times \left({\displaystyle \frac{{\mathit{\mu }}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}}{\mathrm{0.151}}}+\mathrm{0.079}\times {\displaystyle \frac{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}{{\mathrm{N}}_{\mathrm{tot}}}}\right)\times B\end{array}$$

The anaerobic fraction anvf is described following Eq. (13):

$$\begin{array}{}\text{(13)}& \text{anvf}=\mathrm{0.85}\times \left(\mathrm{1}-{\displaystyle \frac{p_{\text{soil}{\mathrm{O}}_{\mathrm{2}}}}{p_{\text{air}{\mathrm{O}}_{\mathrm{2}}}}}\right),\end{array}$$

with $p_{\text{air}{\mathrm{O}}_{\mathrm{2}}}$ and $p_{\text{soil}{\mathrm{O}}_{\mathrm{2}}}$ being the partial pressure in the air and in the soil, respectively. $p_{\text{soil}{\mathrm{O}}_{\mathrm{2}}}$ is calculated following Eq. (14).

$$\begin{array}{}\text{(14)}& {\displaystyle \frac{\partial p_{\text{soil}{\mathrm{O}}_{\mathrm{2}}}}{\partial t}}=p_{\text{soil}{\mathrm{O}}_{\mathrm{2}}}-p_{{\mathrm{O}}_{\mathrm{2}}\text{resp}}\times k\times \text{SOC}\times f_{\mathrm{Cu}}\end{array}$$

Here SOC is the soil organic carbon stock (g C m⁻²), k the decomposition rate, $p_{{\mathrm{O}}_{\mathrm{2}}\text{resp}}$ the O₂ partial pressure related to the respiration and f_Cu the effect of Cu on CO₂ emissions as defined in Eq. (15), following (Sereni et al., 2021; Eq. 5)

$$\begin{array}{}\text{(15)}& f_{\mathrm{Cu}{\mathrm{CO}}_{\mathrm{2}}}=\exp (-\mathrm{0.1}-\mathrm{0.1}\times \log (\mathrm{Cu})+\mathrm{0.12}\times \text{pH}).\end{array}$$

The relative growth rate of NO₃-N, NO_x-N and N₂O-N denitrifiers is described, respectively, by ${\mathit{\mu }}_{{\mathrm{NO}}_{\mathrm{3}}}$, ${\mathit{\mu }}_{{\mathrm{NO}}_{\mathit{\text{x}}}}$ and ${\mathit{\mu }}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}$ following Eqs. (16)–(18).

$$\begin{array}{}\text{(16)}& {\displaystyle }& {\displaystyle }{\mathit{\mu }}_{{\mathrm{NO}}_{\mathrm{3}}}={\displaystyle \frac{\mathrm{0.67}\times f_{\text{denit}}\left(\text{temp}\right)\times f_{{\text{denit}}_{{\mathrm{NO}}_{\mathrm{3}}}}\left(\text{pH}\right)\times {\mathrm{NO}}_{\mathrm{3}}}{{\mathrm{NO}}_{\mathrm{3}}+\mathrm{166}}}\text{(17)}& {\displaystyle }& {\displaystyle }{\mathit{\mu }}_{{\mathrm{NO}}_{\mathit{\text{x}}}}={\displaystyle \frac{\mathrm{0.34}\times f_{\text{denit}}\left(\text{temp}\right)\times f_{{\text{denit}}_{{\mathrm{NO}}_{\mathit{\text{x}}}}}\left(\text{pH}\right)\times {\mathrm{NO}}_{\mathit{\text{x}}}}{{\mathrm{NO}}_x+\mathrm{166}}}\text{(18)}& {\displaystyle }& {\displaystyle }{\mathit{\mu }}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}={\displaystyle \frac{\mathrm{0.34}\times f_{\text{denit}}\left(\text{temp}\right)\times f_{{\text{denit}}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{0}}}\left(\text{pH}\right)\times {\mathrm{N}}_{\mathrm{2}}\mathrm{O}}{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}+\mathrm{166}}}\end{array}$$

Here f_denit(temp), $f_{{\text{denit}}_{{\mathrm{NO}}_{\mathrm{3}}}}\left(\text{pH}\right)$, $f_{{\text{denit}}_{{\mathrm{NO}}_{\mathit{\text{x}}}}}\left(\text{pH}\right)$, and $f_{{\text{denit}}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}}\left(\text{pH}\right)$ are rate modifiers depending on air temperature and soil pH described in Eqs. (19) to (22).

$$\begin{array}{}\text{(19)}& {\displaystyle }& {\displaystyle }{f}_{\text{denit}}\left(\text{temp}\right)={\mathrm{2}}^{(\text{temp}-\mathrm{22.5})/\mathrm{10}}\text{(20)}& {\displaystyle }& {\displaystyle }{f}_{{\text{denit}}_{{\mathrm{NO}}_{\mathrm{3}}}}\left(\text{pH}\right)=\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{\mathrm{1}+e^{(\mathrm{4.25}\times \text{pH})/\mathrm{0.5}}}}\text{(21)}& {\displaystyle }& {\displaystyle }{f}_{{\text{denit}}_{{\mathrm{NO}}_{\mathit{\text{x}}}}}\left(\text{pH}\right)=\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{\mathrm{1}+e^{(\mathrm{5.25}\times \text{pH})}}}\text{(22)}& {\displaystyle }& {\displaystyle }{f}_{{\text{denit}}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}}\left(\text{pH}\right)=\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{\mathrm{1}+e^{(\mathrm{6.25}\times \text{pH})/\mathrm{1.5}}}}\end{array}$$

The denitrifier biomass dynamic B (kg m⁻²) is described following Eq. (23).

$$\begin{array}{}\text{(23)}& {\displaystyle \frac{\partial B}{\partial t}}=\left(\text{anvf}\times ({\mathit{\mu }}_{{\mathrm{NO}}_{\mathrm{3}}}+{\mathit{\mu }}_{{\mathrm{NO}}_{\mathit{\text{x}}}}+{\mathit{\mu }}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}})-\mathrm{3.82}\times {\mathrm{10}}^{-\mathrm{3}}\right)\times B\end{array}$$

Finally, all the gaseous forms of mineral N are emitted into the atmosphere. It is important to note that we did not directly use the DNDC model but a simplified version adapted by Zaehle and Friend (2010). The original code was in Fortran, and we translated it into R to facilitate its manipulation. The time step of the model was 30 min, and most of the parameters were kept to the original values of Li et al. (1992, 2000), except k_nit, which was modified to 0.1743 instead of 0.2 to better fit the data from the control. Furthermore, the amounts of NH₄-N fixed to the clay were reduced to 0 as the bioassay was performed in excess of NH₄-N (see 2.2.0).

Table 1 N species measured in the soils at the end of initial incubation period further used to initialise the DNDC model, mean modelled NO₃-N stocks and mean emissions of NH₄-N, N₂-N, N₂O-N and NO_x-N modelled without Cu. Notation is as follows: 90 is 90 % WHC, 60 is 60 % WHC, DO is Dry-Only and DR is Dry-Rewet treatment during initial incubation. A, B and C are replicates.

n/a: not applicable.

Figure 2 Fitted functions of potential nitrifying activities (PNA) against total soil copper concentrations [Cu] for each initial moisture incubation treatment. Points are the measured nitrate production and lines the fitted quadratic function with their 95 % confidence interval. (a) Constant moisture treatments: the green circle is for 30 % WHC, the red square for 60 % WHC and the purple diamond for 90 % WHC. The black line is the common fitting function used for 60 % and 90 % WHC moisture treatments. (b) Variable initial moisture treatments: the brown star is for Dry-Rewet (DR) and the yellow triangle for Dry-Only (DO).

We used measures of N species at the end of the initial incubation period as initial values of N species for DNDC (Table 1a and Fig. 2). To estimate the anaerobic volume fraction during the 3 d bioassay, we used a C mineralisation rate k (Eq. 14) determined on the basis of measurements performed on the same soil (Annabi et al., 2007) and ran DNDC for a 45 d equilibrium period. We then extracted the initial anaerobic volume fraction and partial O₂ pressure.

2.4 Statistical analysis

The dose–response curves of PNA during the bioassay to Cu gradient were plotted and tested with linear, quadratic or cubic functions as fitting models. Our aim was to find, if possible, a similar modelling fit function for all initial moisture incubation treatments. Thus, for each moisture treatment, the two best functions of fit were selected through Akaike information criterion (AIC) and R² criteria and compared with ANOVA. After selection of a common type of functions, the permutability of the different function parameters was tested with the Chow test (gap v.1.2.2 package, which tested regression 1 on the basis of sample 2 and vice versa). If the p.v. (p value) exceeds its critical values, regressions cannot be considered equal (Zhao, 2007).

To estimate the effect of [Cu] and soil moisture on the different variables measured, we used the nonparametric Kruskal–Wallis test. The fits between the model and the data of soil nitrate concentration during the bioassays were measured using root-mean-square error (RMSE, Eq. 2):

$$\begin{array}{}\text{(24)}& \text{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{N}}{\sum }_{i=\mathrm{1}}^N(X_i-Y_i{)}^{\mathrm{2}}},\end{array}$$

where i is the number of observations (1 to N), X is the predicted value and Y is the observed value. RMSE was decomposed in standard bias (Eq. 25), non-unity slope (Eq. 26) and lack of correlation (Eq. 27) components following Gauch et al. (2003), with $\stackrel{\mathrm{‾}}{X}$ and $\stackrel{\mathrm{‾}}{Y}$ the mean modelled and observed values, b the slope of the least-square regression of Y on X and r² the square of the correlation.

$$\begin{array}{}\text{(25)}& {\displaystyle }& {\displaystyle }\text{SB}={\left(\stackrel{\mathrm{‾}}{X}-\stackrel{\mathrm{‾}}{Y}\right)}^{\mathrm{2}}\text{(26)}& {\displaystyle }& {\displaystyle }\text{NU}=(\mathrm{1}-b{)}^{\mathrm{2}}\times \sum {\displaystyle \frac{x_n^{\mathrm{2}}}{N}}\text{(27)}& {\displaystyle }& {\displaystyle }\text{LC}=(\mathrm{1}-r^{\mathrm{2}})\times \sum {\displaystyle \frac{y_n^{\mathrm{2}}}{N}}\end{array}$$

All analyses were done with R 3.2.3 (R Core Team, 2015).

3 Results

3.1 Effect of Cu on potential nitrification activity (PNA): statistical model selection

The soil N species measured at the end of the soil initial incubations in each soil moisture treatment were used to initialise the DNDC model (Table 1). Two anomalous points leading to anomalous calculated NO₂-N values were excluded from the experimental results because of technical problems during measurements (the C replicates in the DR and DO cases).

The bioassay experiments performed at the end of the soil initial incubations allowed us to determine the rate of nitrate production as a function of soil [Cu] for each soil moisture level (Fig. 1). In all cases, the PNA values were found to decrease with the increase in soil [Cu] but at different rates depending on the moisture treatment. Based on AIC values (Table S2 in the Supplement), we first selected one model per moisture incubation that better fit the data. For 30 % and 60 % of WHC, a quadratic model was found to provide the better compromise between the number of parameters and the prediction capacity. For 90 % WHC, no significant difference was found between the cubic and the quadratic models (ANOVA, p.v = 0.07). For DR, no significant difference was found between linear and quadratic models (Table S2a and b), whereas for DO the cubic model provided a substantially better fit than the quadratic model (AIC and adj. R² score, Table S2). Finally, we found that the quadratic model fitted all the sets of data correctly, allowing homogeneity across the initial moisture incubation treatments (Fig. 2b). The quadratic function was thus chosen to quantify the Cu effect on PNA including the DO treatment.

The parameters of the five quadratic functions (one for each moisture treatment) were found to be different from each other, except for 60 % and 90 % WHC (p.v. = 0.001, Chow test). A single function was thus used to adjust PNA to soil [Cu] curves at 60 % and 90 % WHC but with different intercepts for these two WHC treatments (Table S3 in the Supplement and Fig. 2).

The final four quadratic equations are as follow: Eq. (28) for 30 % WHC, Eq. (29) for 60 % and 90 % WHC, Eq. (30) for DR, and Eq. (31) for DO (Fig. 2).

$$\begin{array}{}\text{(28)}& {\displaystyle }& {\displaystyle }{F}_{\mathrm{Cu}\mathrm{30}}=\mathrm{0.782}-\mathrm{0.000451}\times \mathrm{Cu}+\mathrm{9.49}\times {\mathrm{10}}^{-\mathrm{8}}\times {\mathrm{Cu}}^{\mathrm{2}}\text{(29)}& {\displaystyle }& {\displaystyle }{F}_{\mathrm{Cu}\mathrm{60}/\mathrm{90}}=b-\mathrm{0.000342}\times \mathrm{Cu}+\mathrm{4.30}\times {\mathrm{10}}^{-\mathrm{8}}\times {\mathrm{Cu}}^{\mathrm{2}}\end{array}$$

with b = 0.795 for 60 % WHC and b = 0.796 for 90 % WHC

$$\begin{array}{}\text{(30)}& {\displaystyle }& {\displaystyle }{F}_{\mathrm{Cu}\text{DR}}=\mathrm{0.552}-\mathrm{0.000164}\times \mathrm{Cu}+\mathrm{6.09}\times {\mathrm{10}}^{-\mathrm{8}}\times {\mathrm{Cu}}^{\mathrm{2}}\text{(31)}& {\displaystyle }& {\displaystyle }{F}_{\mathrm{Cu}\text{DO}}=\mathrm{0.625}-\mathrm{0.000192}\times \mathrm{Cu}+\mathrm{2.82}\times {\mathrm{10}}^{-\mathrm{8}}\times {\mathrm{Cu}}^{\mathrm{2}}\end{array}$$

According to the fitted equations, the decrease in nitrate production rates as a function of soil [Cu] depended on initial incubation treatment. Decreases were found to be steeper following 30 % WHC > 60 %–90 % WHC > DO > DR.

These four equations were then added in the DNDC model, allowing us to adjust Eq. (2), which regulates the nitrate production to soil Cu contents.

3.2 Modelling soil nitrate concentrations in Cu-contaminated treatments using the DNDC-Cu model

3.2.1 Setup of the DNDC-Cu model

The DNDC model was originally constructed to model both C and N soil cycles. The relative proportion of nitrification and denitrification processes thus depends on soil aerobic fraction determined by both soil C respiration and soil moisture (Eqs. 13 and 14). Before any addition of Cu function in DNDC, we estimated this soil aerobic fraction using C mineralisation. Previous data from 366 d incubations made on the same uncontaminated soil (Annabi et al., 2007) were first used to fit a C mineralisation coefficient rate, k. The resulting k coefficient (k = 1.234 × 10⁻⁴ $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}\mathrm{30}{\min }^{-\mathrm{1}}$) was introduced in the DNDC model and forced to equilibrium (45 d) without soil Cu contamination effect. This provided a basal aerobic volume fraction for each soil moisture through Eq. (13), corresponding to 3.52 × 10⁻³ at 30 %, 6.167 × 10⁻³ at 60 % (and DR–DO to which bioassays were performed at 60 % WHC) and 2.705 × 10⁻² at 90 % of the WHC. The partial O₂ pressure was calculated as 211.4 hPa at 30 % WHC; 210.7 hPa at 60 % WHC, DR and DO; and 205.4 hPa at 90 % WHC. These values were used to initiate the DNDC-Cu version. We then ran the DNDC-Cu version for a 3 d simulation. The constant rate of C mineralisation, k, was adjusted to take into account the Cu contents with Eq. (14) while Eqs. (28)–(31) were used to adjust NO₃-N production rate (Fig. 1) to Cu.

Figure 3 Comparison between modelled and measured soil nitrate concentration incubated in different moisture levels with 1:1 line: (a) the three initial incubations under constant moisture and (b) the two initial incubations under variable moisture Dry-Rewet (DR) and Dry-Only (DO) treatments. For 30 % WHC, model = 1.84 ⋅ measure and R² = 0.93; for 60 % WHC model = 0.93 ⋅ measure. R² = 0.99; for 90 % WHC model = 0.90 ⋅ measure. R² = 0.99; for Dry-Rewet (DR) model = 0.96 ⋅ measure. R² = 0.98; for Dry-Only (DO) model = 0.95 ⋅ measure. R² = 0.99.

3.2.2 DNDC-Cu model validation

Our DNDC-Cu model has been evaluated by comparing experimental data of soil nitrate concentration measured after 1 and 3 d of the bioassay incubation with the model outputs. A good fit was provided for 60 % and 90 % of WHC in the range of the DNDC calibration compared to 30 % WHC where the nitrate production is largely underestimated (by a factor of 2 after 3 d of incubation, Fig. 3a). The regression slopes between modelled and measured soil nitrate concentration for 60 % and 90 % WHC were, respectively, 0.94 ± 0.01 and 0.91 ± 0.01 (R² = 0.99 in both cases, Fig. 3a), whereas for 30 % WHC the regression slope was 1.21 ± 0.08 (R² = 0.92) (Fig. 3a). For DR, the soil nitrate stocks were either overestimated (at 762 mg Cu kg⁻¹ of soil) or underestimated (at 2012 mg Cu kg⁻¹ of soil, Fig. 3b), but overall modelling adequately fit the data with a regression slope at 0.95 ± 0.02 and R² = 0.99. For DO, the regression slope between modelled and measured soil nitrate stocks was 0.95 ± 0.02 too. The Fig. S1 in the Supplement shows the improvement of the DNDC-Cu version to model NO₃-N soil concentration for contaminated soils, with the differences between modelled and measured [NO₃-N] using the default DNDC version compared to our DNDC-Cu version for each [Cu].

Considering all the moisture treatments, RMSE was about 57.3 as a mean (46.4 g NO₃-N m⁻² standard error) for a mean soil nitrate measured at 390 g NO₃-N m⁻² (69 g NO₃-N m⁻² standard error) after 3 d of incubation. However, for the 30 % WHC, RMSE was 139.9, thus 3.7 times more than for the other treatments (Fig. S2 in the Supplement). Despite the reduction in nitrate production rate from 0.20 to 0.18 g N h⁻¹ (see Material and methods), soil nitrate stock was still slightly overestimated in the 90 % WHC as shown by the largest lack of correlation in this case compared to the 60 % WHC treatment (Figs. 3a and S2). Lack of correlation was reduced for all tested moisture treatments (mean $\sqrt{\text{LC}}$ = 23.0, standard error = 5.4, which is roughly $\mathrm{1}/\mathrm{20}$ of the produced nitrate in 3 d in the uncontaminated treatment). Results showed that our DNDC-Cu version was able to reproduce the variability observed in Cu-contaminated soils except for the 30 % WHC treatment where soil nitrate stocks were largely underestimated. The following results thus focused on the use of DNDC-Cu for DR, DO, and 60 % and 90 % WHC treatments to predict soil N emissions.

3.3 Use of DNDC-Cu to predict N fluxes in contaminated soils

3.3.1 Effect of soil [Cu] on soil N stocks

The soil Cu function we included in the DNDC-Cu model specifically modified the default nitrification equation in response to pH, soil moisture and O₂ availability (Eq. 2). In the presence of low [Cu] (12–512 mg Cu kg⁻¹ of soil), the predicted NO₃-N soil stocks were found to be equivalent between 60 % WHC and DO and, to a less extent, DR treatments (Fig. S3 in the Supplement). When soil [Cu] increased, soil [NO₃-N] decreased but with different rates depending on the moisture of initial incubations (Eqs. 28–31). The evolutions of concentrations in soils and emissions fluxes of each species in response to [Cu] gradient were also found highly different depending on the species and on the moisture of initial incubations. However, the relative evolution in terms of both soil concentration and emissions fluxes was identical for each species and each initial incubation treatment and is represented in Table 2. The largest variations were modelled for N₂O-N decrease (around −63 % for the constant moisture treatments and −54 % for the DR at 2012 mg Cu kg⁻¹ of soil) while the smallest variations were modelled for NH₄-N increase (8 %–10 % for the 60 % and 90 % WHC against 5 %–7 % for the DR and DO initially incubated soils at 2012 mg Cu kg⁻¹ of soil). Due to the different evolutions with Cu gradient, concentrations or intensities of fluxes for a given species may reverse between two moisture treatments with an increase in soil [Cu].

For instance, up to 548 mg Cu kg⁻¹ of soil, we modelled the lowest NO₃-N stocks in DR incubated soils. Above this level, NO₃-N soil stocks were the smallest for the 60 % WHC treatment as a result of the sharpest decrease in NO₃-N production due to soil [Cu]. NO₃-N soil stock for initial incubation at 90 % WHC was the highest for soil [Cu] below 1432 mg Cu kg⁻¹ of soil. Between 1432 and 2000 mg Cu kg⁻¹ of soil, NO₃-N soil stocks were similar for 90 % WHC, DR and DO (Fig. S3).

Table 2 Percentage of variation in soil NO₃-N stocks; soil ${\mathrm{NO}}_{\mathrm{3}}\text{-}\mathrm{N}/{\mathrm{NH}}_{\mathrm{4}}\text{-}\mathrm{N}$ stocks; and NH₄-N, N₂-N, NO_x-N and N₂O-N emissions in response to soil [Cu] in the various initial incubation treatments for a 3 d modellisation.

In the absence of Cu, ${\mathrm{NO}}_{\mathrm{3}}\text{-}\mathrm{N}/{\mathrm{NH}}_{\mathrm{4}}\text{-}\mathrm{N}$ ratios were similar among soil moisture treatments. However, the variations in NH₄-N and NO₃-N stocks in response to Cu gradient were different across soil moisture levels. Indeed, the increase in soil [Cu] resulted in a decrease in nitrification rate and thus in an increase in soil NH₄-N stocks (Fig. S4 in the Supplement). The ${\mathrm{NO}}_{\mathrm{3}}\text{-}\mathrm{N}/{\mathrm{NH}}_{\mathrm{4}}\text{-}\mathrm{N}$ stock ratios decreased faster for 60 %–90 % WHC than for DR and DO with an increase in soil [Cu] (Fig. S5 in the Supplement, Table 2).

The decrease in soil NO₃-N stocks at high [Cu] induced a decrease in the modelled growth of denitrifying bacteria that is directly related to [NO₃-N] (Eq. 13). Consequently, the modelled denitrifying bacterial pool was reduced when soil [Cu] increases (Fig. 4). Whatever the soil [Cu], denitrification was modelled as being larger by roughly a factor of 2 in the soils incubated at 90 % WHC compared to the other treatment as this moist treatment is defined as the perfect condition for denitrifying bacteria in the DNDC model (Li et al., 1992). Soils incubated at 60 % WHC were modelled with the lowest denitrifying bacterial pool. No difference between the DR and DO soils was found due to uncertainties in the modelled denitrifying bacterial pool, which resulted from the different concentrations in N species used to initialise DNDC-Cu (Table 1). The soil N₂O-N stocks and dissolved NO_x-N being directly related to denitrifying bacteria, they followed trends similar to those of soil NO₃-N stocks with a global decrease in soil stocks with an increase in soil [Cu] (table 2) and larger stocks at the wetter treatment.

Figure 4 Modelled soil denitrifying bacterial pool after 3 d (kg m⁻³ soil) for the four moisture treatments. The purple diamond is for 90 % WHC, red square for 60 % WHC, brown star for Dry-Rewet (DR) and yellow triangle for Dry-Only (DO). Red, brown and yellow curves are superposed. Pools were modelled for 12, 62, 112, 262, 512, 762, 1012 and 2012 mg Cu kg⁻¹ of soil as represented by stars. Quadratic fits were used for representation.

Figure 5 Modelled N emission fluxes at 3 d in $\mathrm{g}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}\mathrm{30}{\min }^{-\mathrm{1}}$ under the different moisture treatments. (a) NH₄-N emission fluxes. (b) NO_x-N emission flues, (c) N₂O-N emission fluxes and (d) N₂-N emission fluxes. The purple diamond is for 90 % WHC, red square for 60 % WHC, brown star for Dry-Rewet (DR) and yellow triangle for Dry-Only (DO). Fluxes were modelled for 12, 62, 112, 262, 512, 762, 1012 and 2012 mg Cu kg⁻¹ of soil as represented by stars. Quadratic fits were used for representation.

3.3.2 Estimation of soil N emissions under various moisture levels

Large differences are predicted in the NH₄-N, NO_x-N and N₂O-N fluxes between the 90 % WHC soil and the three other soil moisture treatments (Fig. 5). Due to the different evolutions of fluxes in response to Cu, NH₄-N fluxes were modelled as being the smallest for the DR soils compared to the 60 % WHC incubated for soil Cu below 1774 mg Cu kg⁻¹ of soil and higher above 1774 mg Cu kg⁻¹ of soil (Fig. 5a). The emissions of NH₄-N in the DO treatment were predicted to be higher than those of the DR treatment for soil Cu higher than 1290 mg Cu kg⁻¹ of soil and the smallest below 1290 mg Cu kg⁻¹ (Fig. 5a). In the studied range of added Cu, NO_x-N fluxes predicted by the model are largest from 60 % WHC to DO, DR and 90 % WHC (Fig. 5b) for a moderate Cu input (∼ below 1380 mg Cu kg⁻¹ of soil). The decrease in NO_x-N emission with the increase in soil [Cu] was however steeper for soils incubated at 60 % WHC (Tables 2a and b). Hence, at 2012 mg Cu kg⁻¹ of soil NO_x-N fluxes in soil incubated at 60 % WHC were similar to those in the soils incubated under drought treatment (Fig. 5b). The smallest fluxes of N₂O-N were predicted for the wetter treatment despite higher modelled N₂O-N stocks at 90 % WHC whatever the [Cu] (Table 2a and Fig. 5c). The N₂O-N emissions fluxes in the presence of Cu were predicted to be 4 times smaller in the 90 % WHC treatment compared to the others. N₂O-N fluxes had similar trends to NO_x-N for moderate Cu inputs, but fluxes were still the largest from 60 % WHC to DO, DR and 90 % WHC (Fig. 5c), and N₂-N emissions were larger in the wettest treatment (Fig. 5d). The ratio of emitted N₂O-N per denitrification product (i.e. ${\mathrm{N}}_{\mathrm{2}}\mathrm{O}\text{-}\mathrm{N}/{\mathrm{N}}_{\mathrm{2}}\mathrm{O}\text{-}\mathrm{N}+{\mathrm{N}}_{\mathrm{2}}-\mathrm{N}$) was hence the smallest in the moistest soils, which means that the largest soil N₂O-N stocks in the case of 90 % WHC had more chance to be transformed rather than emitted (Fig. 6).

Figure 6 Proportion of N₂O-N emitted arising from the denitrification calculated as ${\mathrm{N}}_{\mathrm{2}}\mathrm{O}\text{-}\mathrm{N}/{\mathrm{N}}_{\mathrm{2}}\mathrm{O}\text{-}\mathrm{N}+{\mathrm{N}}_{\mathrm{2}}-\mathrm{N}$ modelled fluxes in response to soil Cu concentration for the various moisture treatments. The red square is for 60 % WHC, purple diamond is for 90 % WHC, yellow circle for Dry-Rewet (DR) and brown star for Dry-Only (DO). Red, yellow and brown curves are superposed. Fluxes were modelled for 12, 62, 112, 262, 512, 762, 1012 and 2012 mg Cu kg⁻¹ of soil as represented by stars. Quadratic fits were used for representation.

4 Discussion

4.1 From laboratory experiment to soil N emission modelling

Thanks to our laboratory experiments, we were able to define a function modulating the soil NO₃-N production rates in relation with soil [Cu] and depending on soil moisture. Our results showed that soil nitrate decreases with an increase in soil [Cu]. Initial incubation treatment significantly affects the response of soil nitrate production rate to subsequent Cu stress, with a steeper decrease on the order 30 % WHC > 60 %–90 % WHC > DO > DR for the Cu range studied. The lowest sensitivity of Cu in soils initially incubated with dry–rewet events suggests that it might have selected more resistant communities (Barnard et al., 2013; Gleeson et al., 2008). More complex dose–response functions have been used in Sereni et al. (2022) to assess thresholds and loss of function after such a double stress. These results are in relatively good agreement with those presented here using the quadratic fit, especially for the highest half of [Cu]. However, they also presented a limited increase in nitrification rate for small Cu input that we were not able to emphasise in the present study. In the present article we used simple functions of fit to describe the response of soil nitrate production to Cu gradient after the first moisture stress as they further have to be included in the DNDC model. After implementing these quadratic Cu modulating functions into the DNDC-Cu model, we were however able to reproduce the observed soil nitrate stock, particularly for the soils incubated at 60 % and 90 % WHC. The variability around the mean due to the Cu effect was also reproduced by our DNDC-Cu version at 30 % of WHC despite strong underestimation of mean soil nitrate stocks due to model moisture limit (Li et al., 1992). In the case of the DR and DO incubated soils, the so-called “Cu function” also accounted for the effect of drought stress. In fact, our Cu functions were defined on the basis of soil nitrate production against the whole gradient of Cu, thus also considering the control without Cu. However, the double-stress effect was also well reproduced in nitrate production.

4.2 Expected ecological implications of soil Cu contamination

Based on nitrate production measurements, we modelled a decrease in denitrifying activities with an increase in soil [Cu] as a consequence of the decrease in soil nitrate stocks. However, the experiments performed here did not allow us to determine if the soil Cu contamination rather affects nitrifying bacteria (e.g. decrease in nitrifying activity and in NO₃-N production) or denitrifying bacteria (e.g. increase in denitrifying activities and NO₃-N consumption). The effect of soil contamination on N₂O-N production is debated because (i) microbial species involved are not clearly identified (Wrage-Mönnig et al., 2018), (ii) species richness is not necessarily related to soil functions (Ruyters et al., 2013) and (iii) denitrifying communities could have different sensitivities than nitrification to soil contamination (Hund-Rinke and Simon, 2008; Vásquez-Murrieta et al., 2006). Also, our approach to model N₂O-N, N₂-N and NO_x-N production in the contaminated context could have been more constrained with measurement of denitrification rate to assess the effect of Cu on proportion of production and consumption of NO₃-N.

Based on our simulations, the soil Cu contamination was expected to substantially modify the proportion of available N in soils with the increase in NH₄-N stock at the expense of NO₃-N. NH₄-N accumulation and the large expected decrease in ${\mathrm{NO}}_{\mathrm{3}}\text{-}\mathrm{N}/{\mathrm{NH}}_{\mathrm{4}}\text{-}\mathrm{N}$ ratio in contaminated soils (around 50 % for the 60 % WHC) may lead to a shift in plant community structures with different preferences in N assimilation (Cui and Song, 2007; Peacock et al., 2001). Therefore, Cu stress could not only have implications in microbial community patterns as a stressor but could also induce further shifts due to N species redistributions in soils.

4.3 From N₂O-N, N₂-N and NO_x-N soil stocks to emissions

In the present study, we predicted the highest soil N₂-N, N₂O-N and NO_x-N stocks in the moistest treatments. Indeed these species are produced by the denitrifying bacteria expected to behave optimally at 90 % WHC or after DR cycles (Li et al., 1992; Homyak et al., 2017). However, N₂O-N and NO_x-N emissions were modelled as higher in the driest soils, whereas numerous studies (Dobbie and Smith 2003; Xiong et al. 2007; Manzoni et al. 2012) reported high measured N₂O-N emissions with high moisture. In the present version of DNDC-Cu, the soil N emissions were directly controlled by their diffusion in soil, calculated on the basis of clay and soil moisture content. The diffusion of each species would hence be 11 times smaller under the 90 % WHC (D_s =  0.00357) than under the 60 % WHC treatment (D_s = 0.0306) because the model described the diffusion as a whole and did not separate pores with or without water. Diffusion was hence slower in the water than in the air. Thus, the weighted mean diffusion was lower in the high-moisture treatment. In the absence of Cu soil nitrous stocks being roughly 1.6 times and soil N₂-N stocks 11.1 times larger under 90 % WHC treatment than the other, the emissions of N₂O-N were larger under the driest treatment even if stocks were smaller.

Several studies also reported flushing events with Dry–rewet cycles which would enhance C mineralisation, known as the Birch effect (Birch, 1958; Göransson et al., 2013), hence reducing soil O₂ concentration. Moreover, soil [O₂] is closely related to the pore size distribution, being of major importance in nitrification–denitrification control (Khalil et al., 2004), with a dominating nitrification for aggregates up to 0.25 cm (Kremen et al., 2005). Pore size distribution under dry–rewet events is controlled by cracking, (des)aggregation (Cosentino et al., 2006; Denef et al., 2001) or gas displacement (Kemper et al., 1985) that we were not able to take into account in the present study. In DNDC, the calculation of denitrification rate and diffusion was based on a rough description of the anaerobic zone with approximation of soil pore space distribution (Blagodatsky et al., 2011; Li et al., 2000). The soil pore space distribution approach has been demonstrated to be more generally applicable (Arah and Vinten 1995; Schurgers et al. 2006), whereas soil aggregates have been shown to control the extent of nitrification and denitrification (Kremen et al., 2005; Schlüter et al., 2018). However, if models have been proposed to take O₂ availability at the aggregate size into account in the nitrous oxide production (Kremen et al., 2005; Leffelaar, 1988), they also point out the difficulty in parameterisation that needs a large panel of soil measurements. Moreover, they are rarely transposable at the mesoscale and regional scale due to high spatial variations in soil structure (Butterbach-Bahl et al., 2013). The DNDC-Cu version we used here in particular pointed out the difficulty in dealing with a biogeochemistry model with physical processes, with large discrepancies between modelled soil stocks and emissions. The validation we performed focused on soil nitrate stocks, and a second step to go further on would be the measure of gaseous species to ensure that emissions were also impacted by soil treatment. Moreover, we assumed here that soil [Cu] affected the C mineralisation with a decrease in soil O₂ production leading to an increase in denitrification and N₂O-N and NO_x-N. Nevertheless, the present DNDC-Cu version did not take into account the retroaction between C and N cycles. Further research would thus be required to include Cu contamination in interacting C and N cycles.

4.4 Climate change could substantially modify contaminated soil N emission

It is well known that climate change and rainfall patterns could substantially modify the soil N balance and its GHG emissions (Galloway et al. 2003, 2008; Butterbach-Bahl et al. 2013). Despite limitation in DNDC accuracy for nitrous emissions (Foltz et al., 2019), our results tend to show that increased Cu contamination might also affect soil N emissions, producing the smallest emissions of NO_x-N and N₂O-N. These two gases are of major importance in GHG mitigation with a warming potential per mass 300 and 40 times greater than CO₂, respectively. Agricultural soils being the dominating source of N₂O-N (Beauchamp, 1997; Signor and Cerri, 2013), even a limited decrease in their emissions could have major implications for climate. Based on our modelling, the combined effect of soil moisture and [Cu] was particularly important with larger differences in N₂O-N and NO_x-N emissions between rainfall patterns at high [Cu] (Sect. 3.3.2). Sereni et al. (2022) also showed that soil Cu contamination differently affects soil nitrification depending on primary soil moisture stress. Here we showed that the N₂O-N and NO_x-N emission variations are significantly more sensitive to the combined effect of Cu and precipitation regime than the nitrate stock. Based on these results, Cu inputs to moist soils would lead to larger decreases in soil N₂O-N and NO_x-N emissions compared to drier soils, with the strongest effects likely to occur on soils subjected to abrupt and intense shifts in rainfall patterns.

5 Conclusion

In the present study, we aimed at combining ecotoxicological experiments and biogeochemical modelling focusing on the effect of soil Cu contamination on soil N emission under different soil moisture treatments of constant moisture (30, 60 or 90 % WHC) or a single long drought period (DO) or several dry–rewet (DR) cycles. We showed that the effect of soil Cu contamination was different among moisture treatments and N species. For instance, we modelled that the largest [Cu] (2012 mg Cu kg⁻¹ of soil) provoked a decrease in soil nitrate stocks from −28 % in the DR case to −44 % in the 60 % WHC, whereas N₂O-N emissions were expected to decrease up to 63 % in the 90 % WHC (−62 % in the 60 % WHC case, −54 % in the DO case). However, our results tended to show that the amount of N₂O-N emitted from denitrification would decrease with an increase in soil [Cu] and from 60 % WHC to DR, DO and 90 % WHC, so that less N₂O-N produced would be converted to N₂-N. This result points out two main difficulties in biogeochemical modelling: (i) the difficulty to take into account hydrological dynamics (produced NO₃-N and NH₄-N could be expected to leach) and soil structures at different spatial scales (denitrification is estimated based on rough estimation of anaerobic soil volume which also controlled emissions rates through diffusion processes) and (ii) linking soil function to microbial dynamics, in particular in this case, where only the NO₃-N stock was measured (without dealing between production and consumption for instance). Despite these two main points of uncertainty, the combination of incubations and of modellisation we conducted here emphasises the need to account for soil contamination when dealing with soil GHG emission modelling and climate change, as both contamination and rainfall patterns affect the soil NO_x-N and N₂O-N emissions in a different way.