Study site description
The studied Low Arctic tundra heath is located in a north-south directed valley off the coast of Western Greenland (69°18′40.9″N; 53°30′40.9″W) and described in detail in Rasmussen et al. (2021). During 1991–2018, mean annual air temperature was − 3 ± 1.8°C and mean total precipitation 418 ± 131 mm y− 1 with 42 % falling as snow (Hansen et al. 2006; Zhang et al. 2019). Permafrost underlies most of the area, and on a west-facing slope, a snow fan supplies meltwater downslope during most of the summer. The area has a small external N input of ~ 1 kg N deposition ha− 1 y− 1 and 1–2 kg N fixation ha− 1 y− 1 (Hobara et al. 2006; Rousk et al. 2017).
At the footslope, vegetation is dominated by evergreen and deciduous dwarf shrubs including Betula nana, Empetrum nirgum, Cassiope tetragona, Salix glauca, Vaccinum uliginosum, interlain by mosses such as Tomentypnum nitens, Racomitrium lanuginosum, Sphagnum sp. and lichen. In 2018, a study site was established on this footslope where water drains in the shallow thaw layer during snowmelt and deeper in the later parts of the growing season (Rasmussen et al. 2021). The study site soil temperature, soil moisture, soil water chemistry, plant C:N, and surface gas exchange of N2O was monitored closely throughout 2018 and 2019 (Fig. 1; Rasmussen et al. 2021).
Further, soil was analyzed for grain size distribution, vertical and horizontal hydraulic conductivity to 60 cm depth, C and N content in bulk soil, microbes and roots 0–30 cm, and N2 and N2O development was measured in soil cores sampled in the early and late growing season (Rasmussen et al. 2021). Aboveground vegetation was analyzed for C and N in stem and leaf, split into different plant types (Rasmussen et al. 2021).
Model description and setup
CoupModel is a process-based, numerical ecosystem model that simulates a one dimensional soil-plant-atmosphere system, where biological and abiotic processes are fully coupled via water, heat and mass transfer (Jansson and Karlberg 2011; Jansson 2012). CoupModel has successfully been used to simulate energy and water movement in the Arctic (Hollesen et al. 2011; Rasmussen et al. 2018) and arctic tundra C dynamics (Zhang et al. 2019). Nitrogen turnover, plant-soil interaction and drainage were built in to the predecessor of the CoupModel (Johnsson et al. 1987), and have been applied to simulate grassland (Korsaeth et al. 2003) and forest ecosystems (Christiansen et al. 2006). A trace gas model was later built into the CoupModel and used successfully for agricultural and forest applications as well (Norman et al. 2008; He et al. 2016).
In this study, the model structure was based on a 25-layer soil profile, above which three plant types corresponding to evergreen dwarf shrubs, deciduous dwarf shrubs and moss were simulated. Based on the measurements of soil grain size distribution in the 0–60 cm soil, soil water and temperature, pools of C and N in the 0–30 cm soil and in plant compartments, soil texture and the initial conditions of soil water, temperature, C and N contents of aboveground biomass and root biomass were specified (Table 1). Soil texture for deeper layers (below 60 cm) was based on soil cores from a nearby site (e.g. Zhang et al. 2019). Table S1 shows the model profile with corresponding depths of data collection. Layer-specific initial conditions of soil C and N (Table S2), soil properties related to water retention and hydraulic conductivity (Table S3), C and N contents in aboveground plant compartments (Table S4.1 and S4.2), and soil Total N and C pools and water chemical composition over the summers of 2018 and 2019 (Table S5.1 and S5.2) are also given in the supplementary information.
Table 1
Measured variables specified in the soil profile setup of CoupModel.
Depth
|
Dominating vegetation type (%)
|
Hydraulic conductivity vertically + horizontally ( m s− 1)
|
Soil water chemical composition
(µg L− 1)
|
Soil
Total N
(g g dw− 1)
|
Soil
Total C
(g g dw− 1)
|
N in root biomass
(g g dw− 1)
|
C in root biomass
(g g dw− 1)
|
N in aboveground vegetation
(g g dw− 1)
|
C in aboveground vegetation
(g g dw− 1)
|
N2O flux
(g N m− 2 h− 1)
|
Surface
|
x
|
|
|
|
|
|
|
x
|
X
|
x
|
10 cm
|
|
x
|
x
|
x
|
x
|
x
|
x
|
|
|
|
20 cm
|
|
x
|
x
|
x
|
x
|
x
|
x
|
|
|
|
30 cm
|
|
x
|
|
x
|
x
|
x
|
x
|
|
|
|
40 cm
|
|
x
|
|
|
|
|
|
|
|
|
50 cm
|
|
|
|
|
|
|
|
|
|
|
60 cm
|
|
x
|
|
|
|
|
|
|
|
|
Specific traits of the three plant types, (1) deciduous shrubs, e.g. Betula nana, Salix arctica, (2) evergreen shrubs, e.g. Empetrum nigrum, Cassiope tetragona, and (3) moss, e.g. Racomitrium lanuginosum, Aulacomnium turgidum, sphagnum sp., observed in the plots, were specified in the initial conditions of vegetation setup according to the measured traits, C/N ratio, biomass [g m− 2] and surface coverage [%] of each plant type (Table 2).
Table 2
Specified traits for each plant type (deciduous, evergreen and moss) defined in the setup.
Plant type
|
Max plant height (m)
|
Specific Leaf area (m2 kg− 1)
|
Maximum rooting depth (m)
|
Initial Surface cover (%)
|
Max surface cover (%)
|
Initial C/N Leaf
|
Initial C/N Stem
|
Initial N Leaf
(g m− 2)
|
Initial N Stem
(g m− 2)
|
Initial Leaf C (g m− 2)
|
Initial Stem C
(g m− 2)
|
Deciduous
|
0.15
|
30
|
-0.35
|
20
|
60
|
20
|
50
|
2
|
3
|
40
|
80
|
Evergreen
|
0.1
|
30
|
-0.25
|
30
|
60
|
20
|
75
|
0.7
|
0.8
|
50
|
65
|
Moss
|
0.05
|
30
|
-0.02
|
10
|
40
|
30
|
30
|
3.5
|
0.1
|
50
|
30
|
The following modules in the CoupModel have been switched on to simulate processes that regulate ecosystem abiotic and biological processes and soil gas dynamics (Jansson and Karlberg 2011): soil heat and water transport, snow pack, soil evaporation, plant transpiration, plant growth, soil mineral and organic processes. The details of equations and parameters that regulate all these processes are provided in Table S6. The following description mainly covers the details about how the CoupModel simulates plant growth, soil mineral and organic processes (Jansson and Karlberg 2011).
The plant simulation scheme was set to “Dynamic” so that seasonal development of leaf area index, canopy height, root depth and length were estimated based on plant biomass and allometry. Multiple plants were enabled to grow in the same stand accounting for competition within the plant community. Plant transpiration and soil evaporation were simulated separately using the “explicit big leaf” option, which calculated transpiration using the Penman-Monteith equation. The leaf assimilation of CO2 used the light-use efficiency approach, accounting for responses of temperature, water and C:N ratio of leaf. Plant respiration consisted of both maintenance and growth respiration using the “Q10” type temperature response function. Litter turnover was simulated as steps in which above-ground plant litter moves from first a surface litter pool with no microbial activity into the litter pool, and then to the humus pool. Plant N uptake was simulated based on the demand created by potential growth. Root uptake of N was calculated depending on the N demand, root depth and ratio of nitrate to ammonium. If the actual mineral N uptake was less than the plant demand, organic N uptake (amino acids) can optionally be added to the root uptake by letting plants take up a portion of the litter N pool relative in size to the difference in plant demand and actual mineral N uptake, but with a maximum organic N uptake amount.
For simulating soil decomposition, two C pools (litter and humus) were considered and produced CO2, at a fast and slow rate, respectively. Their decomposition rates were controlled by first order kinetics, temperature and moisture response functions. Microbes were generally simulated implicitly in the soil as opposed to the explicit option, but for nitrification and denitrification rates, the microbe-based option was used. With implicit microbes, there was no pool of microbial N, but they were accounted for in the litter pool of N. The mineralization or immobilization by microbes depended on the respective C/N ratios of the decomposing litter pool. Nitrification rates, understood as transformation of N from the NH4+ pool to the NO3− pool, can be adapted to the ecosystem by adapting the response functions of denitrification to temperature, soil moisture and substrate availability. The simulation of denitrification from NO3− to gaseous N2O and N2 depended on the response functions of denitrification to temperature, soil moisture and substrate availability.
Meteorological drivers
The hourly air temperature, global radiation, longwave radiation, wind speed, relative humidity and precipitation for 2013–2019 were obtained from a nearby weather station (~ 200 m) used as drivers of the simulation. Model spin-up drivers used the meteorological data from 1950–2013 according to Zhang et al. (2019). Table S7 shows the climate parameters used for driving the model setup and the measurement periods and frequencies of data collection.
Calibration and validation
After identifying parameters governing microbial turnover rates and plant organic N uptake, a Monte-Carlo sampling of calibration parameters with 10,000 iterations was done for a multi-criteria based calibration. Based on measured soil temperatures and soil moisture in 10 cm from 2013–2019 with coefficient of determination (R2 > 0.9 for soil temperature and R2 > 0.7 for soil moisture) and Mean Error (ME < 0.5°C and < 5 % vol.) as constrains, an ensemble of 25 runs was accepted as the posterior runs that best represent the studied ecosystem abiotic and biological process. The model performance of the accepted ensemble was further evaluated against the measured pools of root C and N, leaf C and N, stem C and N (scaled to g m− 2), N2O fluxes and soil solution mineral and total dissolved N over time. A list of parameters used for calibration is available in Table S8. Table S9 shows the performance statistics of the selected ensemble after calibration, and Figures S1-4 shows how the model has captured the temporal variations in soil moisture (10, 20, 40 and 60 cm), soil temperature (10, 20, 40 and 60 cm), above- and belowground C and N pools for three species and the soil water chemical composition.
Based on the accepted ensemble of 25 runs as the control runs, we constructed sensitivity experiments testing the impact of lateral water and N input, and the effect of near-surface warming, respectively.
Lateral water input experiments
The control run was made by simulating lateral water input at depth increment 25–35 cm without any N input. Based on this, an ambient scenario was simulated for the years 2012–2019, where the lateral input water had a concentration of NO3− and NH4+ representative of measured concentrations in that depth in the summers 2018 and 2019 (0.03 mg L− 1, Rasmussen et al. 2021).
Additionally, two sensitivity simulations were made for the years 2012–2019, where the NO3− and NH4+ concentrations were 10 × ambient levels and 100 × ambient levels (0.3 and 3 mg L− 1), testing the ecosystem resilience to increased (winter) mineralization and thus N release; the latter considered extreme levels and was rather a test of the model stability than of actual realistic input amounts.
In order to test the effect of near-surface warming compared to increased lateral N input, we made four parallel sensitivity simulations with 0, 0.03, 0.3 and 3 mg N L− 1 in lateral input water, but with surface temperatures increased by 2°C. In total, this yielded one control run and eight sensitivity simulations (Table 3):
Table 3
Overview over experimental simulations with lateral water input in depth 25–35 cm.
Surface temperature
|
0 mg L− 1
|
0.03 mg L− 1
|
0.3 mg L− 1
|
3 mg L− 1
|
Ambient T
|
(Control)
|
(Ambient N)
|
(High N)
|
(Very high N)
|
Near-surface warming (+ 2°C)
|
(Warming)
|
(Warming
+Ambient N)
|
(Warming + High N)
|
(Warming
+Very High N)
|
The input followed the thawing of the soil so that no input was made in a frozen soil.
Partitioning of the lateral N input was calculated by comparing the control run to the experimental runs. The N input, AccNinput, is defined as:
$${Acc}_{Ninput}= {N}_{Fix}+{N}_{dep}+{N}_{lateral}$$
1
where NFix is the amount of N fixated by vegetation [g N m− 2], Ndep is the amount of N deposited on leaves [g N m− 2] and Nlateral is the amount of lateral N input on the frozen surface [g N m− 2].
As such, internal fluxes of N, such as litter N input to the soil, is not counted as an input to the ecosystem.
Similarly, the total N output, AccNoutput, is defined as:
$${Acc}_{Noutput}= {N}_{Emission}+{N}_{Leaching}$$
2
where NEmission is the amount of N lost as gaseous N from the system [g N m− 2] and NLeaching is the amount of N lost as (horizontal) drainage out of the profile.
Internal fluxes of N, such as plant N uptake from the soil, is thus not counted as an output from the ecosystem.
For each experimental run, the total lateral N gained by the system was calculated as the difference in in- and output of N in 2018 subtracted the control run difference of N in- and output in 2018:
$$N gain by {ecosystem}_{ambient}={\left(Ac{c}_{input}-Ac{c}_{output}\right)}_{ambient}-{\left(Ac{c}_{input}-Ac{c}_{output}\right)}_{control}$$
3
The difference for each possible N pathway was then calculated as the difference between experimental and control runs, and the difference was related to the total amount of lateral N gained by the system (here the plant pool):
$$\%lateral{N}_{plan{t}_{ambient}}=\frac{PlantNUptak{e}_{ambient}-PlantNUptak{e}_{Control}}{N gain by ecostyste{m}_{ambient}}*100$$
4
By calculating the amount of N in the different pools this way, the fluxes may also include N from other sources than lateral input water, since e.g. plant N uptake also comes from internal cycling of N and not only external N input. Thus, numbers above 100 % or negative numbers indicate a shift in the balance of fluxes indirectly impacted by the lateral N input, which is the only difference between the control and experimental runs.