2.1 SOLVEG

The one-dimensional multi-layer model SOLVEG consists of four sub-models: atmosphere, soil, vegetation, and radiation within the vegetation canopy as shown in Fig. S1 in the Supplement. The general description is available in Katata (2009), Katata and Ota (2017), Nagai (2004), and Ota et al. (2013). Details of the processes of snow accumulation and melting, freeze–thaw in soil, and grassland vegetation growth and development are described in the Supplement.

In the atmosphere sub-model, one-dimensional diffusion equations are solved between atmospheric layers for horizontal wind speeds, potential temperature, specific humidity, liquid water content of fog, turbulent kinetic energy and length scale (Katata, 2009), and gas and aerosol concentrations (Katata and Ota, 2017). In upper boundary conditions, the variables of horizontal wind speeds, potential temperature, and specific humidity (and liquid water content of the fog, gas and aerosol concentrations, if available) are typically obtained from hourly or half-hourly observational data. For further explanation, see Sect. 2.3. Bulk transfer equations are applied at the lowest layer using the soil surface temperature and specific humidity calculated in the soil sub-model. In the soil sub-model, the soil temperature, volumetric soil water content, and specific humidity in the soil pores are predicted based on heat conduction, mass balance in liquid water, and water vapor diffusion equations, respectively (Katata, 2009). Root water uptake is calculated from the transpiration rate in the vegetation sub-model. For CO₂ concentration in soil, mass conservation equations for liquid and gas phases are solved (Nagai, 2004). Organic matter dynamics are also considered (Ota et al., 2013), such as microbial decomposition and dissolved organic carbon (DOC) leaching in the aboveground litter layer, belowground input of carbon from roots (root litter), and soil organic carbon (SOC) turnover and DOC transport along water flows throughout the soil profile for three SOC pools (active, slow, and passive) with different turnover times.

In the vegetation sub-model, profiles of the leaf temperature, leaf surface water, and the vertical liquid water flux are predicted (Nagai, 2004). The heat budget equation at the leaf surface is solved to predict the leaf temperature using key variables from the atmosphere sub-model combined with the radiation scheme. At the upper boundary of the sub-model, the given precipitation intensity is used for calculating vertical liquid water flux within the canopy based on the surface water budget equation. The CO₂ assimilation rate due to photosynthesis is predicted using Farquhar's formulations (Farquhar et al., 1980) and stomatal resistance. In the radiation sub-model, direct and diffuse downward and upward fluxes of solar and longwave radiation are calculated to obtain the radiation energy input at the canopy layers. Fractions of sunlit and shaded leaves at each canopy layer are computed for the stomatal resistance and energy budget calculations.

A multi-layer snow module is mainly developed based on the Community Land Model (CLM; Oleson et al., 2010) and SNTHERM (Jordan, 1991), while the model is unique in including the gravitational and capillary liquid water flows in the unsaturated snow layer based on van Genuchten's concept of water flow in the unsaturated zone (see Hirashima et al., 2010). In the soil module, freeze–thaw processes in soil based on the freezing-point depression equation (Zhang et al., 2007) are considered in heat conduction and liquid water flow equations.

To simulate the winter-related processes for grassland phenology such as leaf development and senescence due to cold stresses, the relevant scheme in the grass growth model named BASic GRAssland model (BASGRA; Höglind et al., 2016) is coupled with the vegetation sub-model of SOLVEG to simulate vegetation growth. The three main features that characterize plant growth in BASGRA are (1) simulation of source–sink relations where the source consists of both current photosynthesis and remobilization of reserves; (2) simulation of leaf area dynamics and tillering for vegetative and generative tillers; and (3) cold hardening and the effect of physical winter stress factors on tiller survival and plant growth. BASGRA has been well validated by using several experimental datasets of harvestable dry matter of perennial rye grass collected in Europe (Schapendonk et al., 1998) and from five locations in Norway, covering a wide range of agroclimatic regions, day lengths, and soil conditions (Höglind et al., 2016). BASGRA consists of the LINGRA grassland model (Van Oijen et al., 2005) with models for cold hardening and soil physical winter processes, while diurnal CO₂ assimilation is calculated as accumulation of the net assimilation for each time step within the vegetation sub-model (Nagai, 2004) instead of the original scheme of photosynthetic processes in BASGRA. When snow covers grass, no photosynthesis is assumed to occur due to low light availability and only soil respiration is considered. BASGRA uses a so-called “big-leaf” approach (Monteith, 1981), thus predicting the total leaf area index (LAI) of the whole grassland vegetation canopies. Since SOLVEG uses a multi-layer structure of canopies, the profile of leaf area density is obtained from simply dividing total LAI by canopy height (h) by assuming vertically uniformity for all canopy layers. Canopy height, which is not simulated in BASGRA, is calculated by the function of LAI with fitting parameters.

The carbon gain from photosynthesis and remobilized reserves is allocated among sinks based on changing sink priorities and strengths. Sink strengths are calculated based on the dynamics of leaves and stems and the acclimation to low temperature. The following five sinks are considered: the processes of cold hardening, replenishment of the reserve pool, leaf growth, stem growth, and root growth. Sink strengths are defined as the rate at which these processes would proceed with no source limitation. The hardening process has top priority, so its demand is met in full if source strength is large enough, irrespective of the four other sinks. Root growth has lowest priority and depends on carbon unused by other sinks. The strength priority between reserves on the one hand, and leaves and stems on the other hand changes with day length. When day lengths are shorter than a cultivar-specific threshold, reserves have higher priority than stems and leaves, with the opposite during the rest of the year. Leaves and stems have equal priority so they receive carbon according to their sink strengths. The removal of tillers and leaves by cutting can be simulated during the growing season, with subsequent regrowth of the sward. The regrowth rate after cutting is calculated at each phenological stage. Natural turnover of leaves and roots is modeled using typical life spans in years (Arora and Boer, 2005), while BASGRA does not simulate the senescence of elongating tillers or roots. The fraction of roots in soil layers and rooting depth are modeled as a function of root biomass (Arora and Boer, 2003), which may be required to be tested at multiple biomes. Daily amounts of the dead root biomass (root litter) are used as inputs to SOC in the soil sub-model of SOLVEG.

2.2 Empirical parameterizations for cold acclimation

Although the relation between the maximum catalytic capacity of Rubisco (V_cmax) and air temperature is quite well established (e.g., Bernacchi et al., 2001; Leuning, 2002; Smith and Dukes, 2013), parameters related to photosynthesis are still uncertain (Kattge and Knorr, 2007) also for low temperature (Höglind et al., 2011). Thus, in the vegetation sub-model, we introduced the empirical factor for cold stress of grasslands, f_cold, to empirically simulate the reduction of photosynthesis under low air temperature as per the following equations (see also Supplement):

where A_n (µmol m⁻² s⁻¹) is the net CO₂ assimilation rate at each canopy layer, which is calculated by subtracting the leaf respiration rate R_d (µmol m⁻² s⁻¹) from the assimilation rate, w_c (µmol m⁻² s⁻¹) is the limitation by efficiency of the photosynthetic enzyme system (Rubisco), w_e (µmol m⁻² s⁻¹) is the limitation by the absorbed photosynthetically active radiation (PAR), w_s (µmol m⁻² s⁻¹) is the limitation by the capacity of leaves to export the products of photosynthesis, T_a (^∘C) is the daily and vertical mean air temperature for all canopy layers, and T_ph (^∘C) is the threshold air temperature above which grasslands are photosynthetically active. Determination of the value of this threshold temperature is important to avoid the overestimation (mainly from fall to winter) of photosynthesis at a low temperature (Höglind et al., 2011). In the original BASGRA, T_ph is set to 1 ^∘C; that is, V_cmax starts decreasing linearly when T_a drops below 1 ^∘C until it becomes zero at −4 ^∘C. However, in the SOLVEG simulation, since the values of T_ph may change depending on environmental conditions, the value of T_ph is calibrated for each site so that the model reproduces the observed CO₂ flux during the extremely warm winter period.

Table 1Characteristics of past or ongoing CO₂ flux observational sites over grassland ecosystems in European mountains. Snow-free CO₂ uptake (bold font) represents the situation of high negative values of CO₂ flux even during the wintertime (typically from December to February).

Download Print Version | Download XLSX

2.3 Study sites and observational data

The model is applied to two sites of managed grassland named Graswang (47.5708^∘ N, 11.0326^∘ E, 864 m a.s.l.) and Fendt (47.8329^∘ N, 11.0607^∘ E, 595 m a.s.l.) belonging to the TERestrial ENvironmental Observatories (TERENO) network in Germany. General information on the climate and management of the sites is available in Table 1. Both sites are located in the Bavarian Alpine Foreland, in the south of Germany and north of the Alps (Mauder et al., 2013; Zeeman et al., 2017, 2019). The grasses are harvested several times during the growing season defined as the period from April to October.

Half-hourly data of precipitation, atmospheric pressure, horizontal wind speed, air temperature, and humidity, and incoming long- and shortwave radiation were used at the top atmospheric layer as a height of 3.5 m. Data of friction velocity (u*), sensible (H) and latent heat (λE), and CO₂ fluxes ($F_{{\mathrm{CO}}_{\mathrm{2}}}$) observed over the grassland based on the open-path eddy covariance method using a three-dimensional sonic anemometer (CSAT3; Campbell Scientific, USA) and an open-path CO₂∕H₂O gas analyzer (LI-7500; Li-Cor, USA) were used for validation of the simulation results. The net radiation (R_net) over the canopies, soil temperature at 0.05 m in depth, and snow depth were also used to evaluate the simulated surface energy and water balances. Details of the site characteristics and micrometeorological observations are described by Zeeman et al. (2017).

Table 2Simulation settings for the modified SOLVEG at the Fendt and Graswang sites. Abbreviations: DM, dry matter; DW: dry weight.

Download Print Version | Download XLSX

2.4 Calibration and validation procedure

Direct comparisons between the results using the original (SOLVEG only) and integrated models (SOLVEG coupled with BASGRA) are difficult because the vegetation dynamics had been prescribed in the original model, requiring time series of total LAI or leaf biomass data, which is used for evaluation in this study. Thus, we simply focus on the calibration of the integrated model only to investigate the impact of wintertime carbon uptake on grassland dynamics. Parameters used for SOLVEG simulations are summarized in Table 2. The study period was approximately 3 years long from 1 December 2011 to 1 November 2014, and included both normal (2011–2012 and 2012–2013) and extremely warm (2013–2014) winters. Typical values of soil hydrological parameters (e.g., saturated hydraulic conductivity) in the study area are given to SOLVEG runs from the past model study (Hingerl et al., 2016). The set of parameters of BASGRA for typical perennial grass species of timothy in the Nordic region (Höglind et al., 2016) is applied. Grass-cutting events are determined from clear reductions in the CO₂ flux, surface albedo, and phenology camera observations according to Zeeman et al. (2017). The threshold temperature for cold stresses (T_ph in Eq. 2) is manually determined in the simulation for each site to obtain the best agreement between simulated and measured CO₂ flux over the canopy during winter. By changing the T_ph value from the range between 1 and 11 ^∘C with an increment of 2 ^∘C (not shown in the figure), we obtained the best results as T_ph=1 ^∘C and 11 ^∘C for Graswang and Fendt, respectively. The calibration results of daily mean surface fluxes (R_net, H, λE, and $F_{{\mathrm{CO}}_{\mathrm{2}}}$) are statistically evaluated using the mean error (ME), the root-mean-square error (RMSE), intercept and slope of linear regression lines, and Pearson's correlation coefficient (R).

Figure 1Time series for (a, c) calculated (lines) and (b, d) observed (open circles) daily mean net radiation (R_net), sensible heat flux (H), and latent heat flux (λE) at (a–b) Fendt and (c–d) Graswang throughout the study period.

Download

2.5 Scenario determination for sensitivity analysis

To investigate the impact of cold acclimation of grassland vegetation on the CO₂ balance and carbon allocation at mountain grassland ecosystems, two scenarios using the SOLVEG model are defined based on the experimental results of Höglind et al. (2011): “active scenario” (T_ph=1 ^∘C) and “dormant scenario” (T_ph=11 ^∘C). The former indicates that photosynthesis is active during most of the wintertime and photosynthesis takes place even at the low temperature of 1 ^∘C. In contrast, the latter represents a situation where grass physiology is more or less shut down and photosynthesis ceases under the conditions of the relatively high temperature of 11 ^∘C to protect from cold stresses. Both scenarios are adopted for both Graswang and Fendt, covering the same period.

Figure 2Time series for calculated (solid lines) and observed (open symbols) (a, c) daily mean soil temperature at a depth of 0.02 m and snow depth and (b, d) CO₂ flux ($F_{{\mathrm{CO}}_{\mathrm{2}}}$) and leaf area index (LAI) at (a–d) Fendt and (e–h) Graswang throughout the study period. Sudden decreases in calculated LAI in (b, d) represent grass-cutting events.

Download

3.1 Model calibration and validation

Figure 1 shows the temporal changes in simulated and observed daily surface heat fluxes over the grassland at Fendt and Graswang throughout the 3-year study period. The model generally reproduced the typical seasonal changes measured at both sites, for example, low values of the Bowen ratio (H ∕ λE) at Fendt during the growing season (from April to October) and negative sensible heat flux (H) at Fendt in December 2013, as suggested by Zeeman et al. (2017).

Figure 3Scatter diagrams of calculated and observed (a, e) daily mean net radiation (R_net), (b, f) sensible (H) and (c, g) latent (λE) heat, and (d, h) CO₂ fluxes ($F_{{\mathrm{CO}}_{\mathrm{2}}}$) at (a–d) Fendt and (e–h) Graswang for the study period. ME (mean error) = calculations − observations; RMSE: root-mean-square error.

Download

Figure 2 illustrates the time series of modeled and observed daily soil temperature and snow depth at the two sites. Observed changes in soil temperature and snow depth (Fig. 2a, c) were also generally reproduced by the model; for example, when the grassland was covered by the snow at Graswang from December 2012 to February 2013, soil temperature at a depth of 0.02 m remained at almost 0 ^∘C in both observations and calculations (Fig. 2c). Sudden increases in soil temperature during periodical snow-free conditions were also reproduced by the model; this was particularly evident at Fendt during the extremely warm winter of 2013–2014 (Fig. 2a).

Figure 4Time series for calculated leaf biomass and snow depth (blue lines) at (a) Fendt and (b) Graswang from 1 September 2012 until 1 June 2014, in active (red lines) and dormant cases (black lines).

Download

Simulated and observed daily CO₂ fluxes ($F_{{\mathrm{CO}}_{\mathrm{2}}}$) over the canopies and simulated LAI at both sites are presented in Fig. 2b. The model simulation match the observed increase of CO₂ uptake due to regrowth of grassland vegetation as well as the sudden decline after harvests (Fig. 2b and d). No drought stress was apparent in the simulations at both sites during the study period (not shown). During the extremely warm winter from December 2013 to February 2014, a net release of CO₂ at Fendt was also well reproduced by the model (Fig. 2b) using the calibrated value of T_ph=1 ^∘C (Table 2). At Graswang, both observed and simulated CO₂ fluxes were very small and close to zero (Fig. 3d), which was achieved in the model by setting a high threshold temperature for cold acclimation (calibrated as T_ph=11 ^∘C; Table 2).

Scatter diagrams and statistical comparisons of daily energy and CO₂ fluxes at the two sites throughout the study period are presented in Fig. 3. At both sites, the slopes of the regression lines were overall close to unity and values of the intercepts were sufficiently small for R_net, H, and λE. High correlations were also observed between measured and simulated CO₂ fluxes at both sites.

3.2 Sensitivity analysis

Figure 4 illustrates temporal changes in simulated snow depth and leaf biomass obtained for the active and dormant scenarios at both sites for the normal winter (2012–2013) and extremely warm winter (2013–2014). It is obvious that leaf biomass during winter is higher in the active scenarios, mostly because of a higher leaf growth at the end of the vegetation period. In addition, some minor leaf growth also occurs in the snow-free winter periods. Nevertheless, leaf biomasses converge during spring and are similar again at the first cutting event in May/June (Fig. 4b).

Figures 5 and 6 depict the cumulative GPP and ecosystem respiration (RE), mean leaf and root biomass, carbon reserve content (total stock of carbon that can be allocated to any of the plant elements such as leaves, stems, and roots), and LAI, simulated for Fendt and Graswang during winter and spring in 2014, respectively. In the following, we focus on Fendt for illustration of the scenario differences (Fig. 5) but would like to emphasize that the responses are similar at both sites. The differences in absolute values, especially a smaller LAI and less biomass are due to the generally cooler conditions at the Graswang site. It should be noted that such differences in environmental conditions between the sites required the calibration of the value of T_ph for each site (Sect. 2.2). Both GPP and RE were higher in the active scenario than in the dormant one as expected by the model construction (Fig. 5a and b); this was particularly apparent as cumulative GPP differed by a factor of 3 or by approximately 100 gC m⁻² yr⁻¹ (Fig. 5a). Changes in leaf biomass and LAI during the subsequent spring in the active scenario were slightly lower than in the dormant scenario, indicating that a small part of the additional GPP has been used for foliage growth (Fig. 5c and f). In contrast, changes in root biomass (below ground) both during winter and spring in the active scenario were approximately 3 times higher than in the dormant scenario (Fig. 5d). Simulated carbon reserve contents in both winter and spring were similar in the two simulation scenarios (Fig. 5e), showing that the carbon fixed by photosynthesis was immediately allocated to the above- or belowground biomass.