Effects of bryophyte and lichen cover on permafrost soil temperature at large scale

Bryophyte and lichen cover on the forest floor at high latitudes exerts an insulating effect on the ground. In this way, the cover decreases mean annual soil temperature and can protect permafrost soil. Climate change, however, may change bryophyte and lichen cover, with effects on the permafrost state and related carbon balance. It is therefore crucial to predict how the bryophyte and lichen cover will react to environmental change at the global scale. To date, current global land surface models contain only empirical representations of the bryophyte and lichen cover, which makes it impractical to predict the 5 future state and function of bryophytes and lichens. For this reason, we integrate a process-based model of bryophyte and lichen growth into the global land surface model JSBACH. The model simulates bryophyte and lichen cover on upland sites, wetlands . :::::::: Wetlands : are not included. We take into account the dynamic nature of the thermal properties of the bryophyte and lichen cover and their relation to environmental factors. Subsequently, we compare simulations with and without bryophyte and lichen cover to quantify the insulating effect of the organisms on the soil. 10 We find an average cooling effect of the bryophyte and lichen cover of 2.7 K on temperature in the topsoil for the study region ::::: region ::::: north :: of ::: 50◦ : N : under current climate. Locally, a cooling of up to 5.7 K may be reached. Moreover, we show that using a simple, empirical representation of the bryophyte and lichen cover without dynamic properties only results in an average cooling of around 0.5 K. This suggests that a) bryophytes and lichens have a significant impact on soil temperature in high-latitude ecosystems and b) a process-based description of their thermal properties is necessary for a realistic representation 15 of the cooling effect. The advanced land surface scheme including a dynamic bryophyte and lichen model will be the basis for an improved future projection of land-atmosphere heat and carbon exchange.

Subsequently, Chadburn et al. (2015b) extend JULES by an empirical scheme which determines moss "health" as a function of climate variables and consequently allows for the computation of a dynamic moss ground cover. In summary, the use of constant thermal properties and empirical relations in current global land surface models makes it difficult to predict the impact of future bryophyte and lichen functions on permafrost ecosystems.
Here we present a process-based model which simulates productivity and dynamic surface coverage of the bryophyte and 5 lichen ground cover at the global scale. The model represents the organisms' water and ice content and thereby accounts for dynamic thermal properties of the cover. The basis for this model is the stand-alone dynamic non-vascular vegetation model LiBry which has been successfully applied to quantify global productivity by bryophytes and lichens (Porada et al., 2013) as well as estimating their contribution to global biogeochemical cycles (Porada et al., 2014). LiBry is fully integrated into the land surface scheme JSBACH of the Max Planck Institute Earth System Model. JSBACH simulates the carbon balance of 10 vascular plants and the soil at the global scale, it includes a water and energy balance and also a representation of permafrost.
Based on the thermal properties of the near-surface vegetation layer we use JSBACH to quantify the effect of bryophytes and lichens on heat transfer between atmosphere and soil and, therefore, on soil temperature. Since JSBACH does not include a scheme for wetland hydrology, the formation of peatlands cannot be simulated. Hence, LiBry in JSBACH mainly represents bryophyte and lichen growth on upland forest floor sites. 15 We compare the soil thermal regime of a JSBACH simulation including a dynamic bryophyte and lichen layer to a simulation where the thermal properties and surface coverage of this layer are set to constant and also to another simulation with the bryophyte and lichen cover switched off. In this way, we can assess quantitatively the impact of the bryophyte and lichen ground cover on soil temperature at high latitudes. In a next step, our process-based approach can be used to predict the role of bryophytes and lichens for high-latitude ecosystems under climatic change.

Model description
For our analysis we use the global land surface model JSBACH 3.0 (Raddatz et al., 2007;Brovkin et al., 2009), which is a part of the Max Planck Institute Earth System Model (MPI-ESM 1.1). JSBACH uses a process-based approach to simulate both physical and biochemical ecosystem functions, such as the exchange of energy and water at the land surface and carbon fluxes 25 between atmosphere, vegetation and soil which are determined by photosynthesis and respiration. The soil part of JSBACH includes a soil carbon model (Goll et al., 2015), a five-layer soil hydrology scheme (Hagemann and Stacke, 2015) and it has been extended by Ekici et al. (2014) to take into account both several dynamic snow layers and the latent heat of fusion associated with freezing and thawing, which is the basis for simulating permafrost extent and active layer thickness.
To represent bryophyte and lichen ground cover, the process-based non-vascular vegetation model LiBry is integrated into 30 JSBACH, which is described in detail in Porada et al. (2013). LiBry is a stand-alone dynamic global vegetation model that uses climate data, such as radiation, temperature and precipitation to predict photosynthesis, respiration and growth. The model combines approaches used in many dynamic vegetation models, such as the Farquhar photosynthesis scheme, with lichen-and bryophyte-specific processes, e.g. poikilohydry.
LiBry uses the Farquhar scheme (Farquhar and von Caemmerer, 1982) to calculate photosynthesis. Thus, both CO 2 and photosynthetically active radiation (PAR) are limiting factors. The availability of PAR on the ground depends on the shading by trees and it is consequently inversely related to the leaf area index (LAI) of the vegetation above the bryophyte and lichen cover. 5 In JSBACH, the fraction of PAR available at the ground can be directly computed from the simulated LAI of the overlying vegetation (Knorr, 2000) and it is then used as an input for LiBry. To estimate the availability of CO 2 for bryophytes and lichens, the model computes diffusion of CO 2 from the atmosphere into the organisms. The CO 2 diffusivity is decreasing with increased water content of the organisms due to narrowing of diffusion pathways and formation of water films (Cowan et al., 1992). In addition to CO 2 and PAR, the rate of modelled photosynthesis depends on surface temperature, which is calculated 10 by JSBACH from the surface energy balance. Furthermore, the photosynthesis rate is related to the level of metabolic activity of bryophytes and lichens, which is controlled by their water content and which ranges from inactive for dry organisms to fully active for water saturated organisms. This is an important function for predicting the response to climate change as precipitation patterns are projected to change (Pachauri et al., 2014).
The water content of LiBry in JSBACH is computed from the balance of water input and loss at a given water storage 15 capacity. For bryophytes and lichens, water input corresponds to rainfall or snowmelt and water loss takes place in form of evaporation. While JSBACH provides fluxes of rainfall and snowmelt, it does not compute the evaporation flux from ground based vegetation. In JSBACH, evapotranspiration is partitioned only into transpiration from vegetation, bare soil evaporation and evaporation from the interception reservoir. The latter includes water evaporating freely from vegetation surfaces in the canopy and on the ground. Since it is impractical to modify this scheme, we approximate evaporation from bryophytes and 20 lichens by evaporation from the interception reservoir. This potentially neglects the morphological control of the organisms on evaporation. However, free evaporation is more suitable than transpiration by vascular plants to describe water loss of bryophytes and lichens (Nash III, 1996;Proctor, 2000). Hence, the water balance of the simulated bryophytes and lichens is coupled to the interception reservoir in JSBACH. This is done by adding the water storage capacity of bryophytes and lichens to the size of the interception reservoir and setting the water saturation of the bryophytes and lichens equal to the saturation of 25 the interception reservoir.
In addition to processes involving liquid water, freezing and melting of water are taken into account in JSBACH. Consequently, water inside the bryophyte and lichen cover may be partly or completely frozen. Since frozen water occupies pore space, the size of the interception reservoir is reduced by the amount of ice inside the bryophyte and lichen cover as long as ice occurs.

30
The net primary productivity (NPP) of bryophytes and lichens is obtained in LiBry by subtracting respiration from gross photosynthesis, where respiration is simulated as a function of temperature and metabolic activity. Turnover of biomass associated with mortality is also considered, and it is modelled as a function of the protein content of biomass. The net growth of bryophytes and lichens is then determined by the balance of NPP and turnover. By multiplying net growth of simulated bryophytes and lichens with their specific area, the growth in surface coverage is derived in the model. Without processes that 35 reduce the bryophyte and lichen cover, the surfage coverage would increase to 100 % in all areas where net growth is larger than zero in steady-state. Hence, a disturbance cycle is included into LiBry which periodically sets back the surface coverage to a small initial value. The surface coverage in steady-state is then obtained by averaging over a whole disturbance cycle.
LiBry in JSBACH is designed to predict the dynamics of the cover in transient scenarios of climate change. Consequently, the steady-state calculation of the surface coverage from the original LiBry is replaced by a dynamic scheme. This also accounts 5 for potential changes in disturbance frequency.
Once a month, the simulated bryophyte and lichen cover is reduced by a certain fraction due to disturbance, such as fire, and the cover is increased by accumulated positive net growth over this period, which is translated into coverage via the specific area of the organisms. Negative net growth is subtracted from the cover. The dynamic scheme also accounts in a simple way for dispersal and establishment of bryophytes and lichens (e.g. Bohn et al. (2011)): The absolute increase in cover per area 10 of ground is both limited by the already existing cover, which generates the new biomass, as well as the free area that is still available for the growth of new cover. The transition of new cover from the existing cover to the free area may take place in form of spores or via vegetative growth (Grime et al., 1990;Rogers, 1990). The efficiency of these processes is summarised in the dynamic scheme by an "expansion efficiency" (see Fig. 1). Equation 1 shows the balance of cover growth and loss where A is the surface coverage of bryophytes and lichens in [m 2 cover m −2 ground], G is net growth in [m 2 new cover m −2 cover], accumulated over a month, and η e is a dimensionless "expansion efficiency" which is set to 0.85. This value is taken from the global, stand-alone version of LiBry, which has also been updated to a dynamic disturbance scheme. It was selected to obtain a realistic global distribution of surface coverage. D is the fraction of cover lost due to disturbance, it is set to 0.00083 per month.
This value corresponds to a fire return interval of 100 years, which is characteristic for the boreal forest (Bonan and Shugart, 20 1989;Beer et al., 2006;Mouillot and Field, 2005). The two minimum statements in Eq. 1 ensure that the cover increment can neither exceed the available area nor lead to negative cover, while the two maximum statements are used to distinguish between positive and negative net growth. The term A(1.0 − A) describes limitation by both existing cover and free area available for growth.
Not only the scheme for computing surface coverage, but also the representation of the organisms' physiological variation 25 had to be modified to integrate LiBry into JSBACH. In the original version of the model, a Monte-Carlo approach is used to sample broad ranges of possible parameter values. From that, artificial "species" are generated that represent the functional diversity of bryophytes and lichens. The "species" perform differently under a given climate. Those which cannot maintain a cover larger than zero cannot "survive" in the model and the remaining ones are used to compute the productivity. This "species"-based approach is conceptually different from the traditional plant functional type (PFT) approach used in many 30 land surface models, and also in JSBACH. To adapt LiBry to the PFT design of JSBACH, the number of artificial "species" is set to one. This "species" is then parameterized to correspond to boreal, ground-based bryophytes and lichens. A is reduced by disturbance and increases due to net growth G. New bryophyte and lichen cover per unit area is limited by the generating area A, the "expansion efficiency" η and the available free area (1 − A).

Bryophytes & Lichens
In general, PFTs in JSBACH are not able to coexist in the same place since they are represented in distinct, non-overlapping tiles. The tiles specify which fraction of a model grid cell is occupied by a certain PFT, thereby representing sub-gridscale heterogeneity. The bryophyte and lichen PFT, however, should be able to grow in combination with various different PFTs.
It is therefore not limited to a specific tile, but instead it is implemented as an additional layer on top of the soil. This layer is available on tiles which contain trees or grasses, while it is excluded from tiles covered by crops or glaciers. The surface 5 coverage of the bryophyte and lichen PFT may differ between the tiles of a given grid cell. This can result from the differential influences of the tiles' vegetation types on bryophyte and lichen growth, for instance due to differing LAI between the tiles.
Defining the bryophyte and lichen PFT as a layer in JSBACH has implications for the representation of the organisms' hydrological properties. The original version of LiBry computes the water content of the organisms through the balance of water input and loss, where the specific storage capacity for water per biomass is sampled for each "species" from a range 10 of possible values from the literature. These values are not related to the geometry of the bryophyte or lichen. In JSBACH, however, the water storage capacity of a layer is determined by multiplying layer thickness by the porosity of the layer. Hence, to be consistent with JSBACH, the water storage capacity of the bryophyte and lichen PFT is computed from the thickness and porosity of the bryophyte and lichen layer. Based on the study by Soudzilovskaia et al. (2013), we set thickness to 4.5 cm, which corresponds to the median of the measured values. The measurements are based on green and undecomposed brown 15 tissue of bryophyte mats. Porosity is not directly measured in Soudzilovskaia et al. (2013), but they provide values of maximum volumetric moisture. Since we are actually interested in the water storage capacity of the bryophyte and lichen layer, we set the "effective" porosity to 80 %, which is at the higher end of the measured values of volumetric moisture.
In JSBACH, two thermal properties of the bryophyte and lichen layer have to be known to derive its influence on soil temperature: Thermal conductivity and heat capacity. Both are strongly dependent on the relative moisture content of the bryophyte and lichen layer. Moreover, they depend on the state of matter of the water in the bryophytes and lichens, which can be liquid or frozen in various relative amounts. In analogy to the soil layers in JSBACH (Ekici et al., 2014), we write the dependence of the thermal conductivity of the bryophyte and lichen layer, κ, on water content as: where κ o is the thermal conductivity of organic matter which is set to 0.25 [W K −1 m −1 ] (Beringer et al., 2001), κ w is the thermal conductivity of liquid water, κ i is the thermal conductivity of ice and κ d is the thermal conductivity of the dry bryophyte and lichen cover, which is set to 0.05 [W K −1 m −1 ] according to (O'Donnell et al., 2009). v w and v i are the volumetric moisture and ice contents of the bryophyte and lichen cover which are calculated by dividing the absolute water or ice content in [m] by 10 he thickness of the cover. K e is the Kersten number which is calculated for the linear regime as: where is the porosity of the bryophyte and lichen cover. The dependence of the heat capacity of the bryophyte and lichen cover C, on water content is written as: 15 where C o is the heat capacity of organic matter which is set to 2.5E6 [J m −3 K −1 ] (Beringer et al., 2001), C w is the heat capacity of liquid water, C i is the heat capacity of ice and ρ w and ρ i are the densities of liquid water and ice. To simulate the influence of the bryophyte and lichen layer on soil temperature, the layer and its thermal properties are included in the vertical heat transfer scheme of JSBACH. The original heat transfer scheme for permafrost soil, which is described in detail in Ekici et al. (2014), already contains an organic layer which represents bryophyte and lichen ground cover. 25 This layer, however, has constant thermal properties since it does not consider the water content of bryophytes and lichens.
It is therefore replaced by the new layer, which explicitly simulates water and ice content of bryophytes and lichens and the associated dynamic thermal conductivity and heat capacity (Fig. 2).
The scheme subdivides the soil column into several layers, with the bryophyte and lichen layer on top. Additionally, a set of snow layers is simulated above the other layers in case snow is present. Thermal conductivity and heat capacity of each layer   are then used to determine the vertical temperature profile by solving the equation for heat conduction for all layers. Thereby, the scheme determines the temperature of a layer in a given grid cell by first calculating the temperature profile for all tiles of the grid cell separately and then averaging the temperatures in each layer weighted by the area fraction of the tiles. The lower boundary condition of the scheme is set by assuming zero heat flux at the bottom of the soil column. The upper boundary condition is surface temperature, which is calculated from the surface energy balance using radiative forcing and the ground 5 heat flux from the previous time step. The scheme also considers the influence of phase change of water on the temperature of a layer, thereby allowing to compute freezing and melting of water in the bryophyte and lichen layer.
Another important difference between the old and the new version of the heat transfer scheme is dynamic surface coverage of bryophytes and lichens. While the surface coverage of the old organic layer is 100 % everywhere, the coverage of the new bryophyte and lichen layer varies between 0 and 100 % between the grid cells of JSBACH (Eq. 1). This means that each 10 tile in a grid cell has a part where the bryophyte and lichen layer is present and another part where the layer is absent. The vertical heat transfer scheme, however, requires a constant number of layers for each tile. Hence, in the new version of the scheme, the temperature profile is calculated twice for each tile, one time for all layers including the bryophyte and lichen layer, and another time leaving out the bryophyte and lichen layer. Subsequently, the average temperature of each layer of the respective tile is obtained by weighting the two profiles with their associated surface coverage. However, the two profiles have 15 a different number of layers, not only due to the bryophyte and lichen layer, but also because the number of snow layers may differ between the parts. Consequently, we fill up the "empty" layers on top by using the surface temperature for the averaging procedure, according to Figure 3.
We implemented two further changes to the heat transfer scheme described in Ekici et al. (2014): We increased the number of soil layers from 5 to 7, thereby extending the layer boundaries at 6.5, 32, 123, 413 and 983 cm depth by boundaries at 23.04 and 53.18 m depth. This ensures that biases due to the zero heat flux lower boundary conditions are negligible and temperature of the lowest layer is not fluctuating after spin-up. Moreover, we complement the snow scheme by dynamic formulations of snow density (Verseghy, 1991) and thermal conductivity (Goodrich, 1982), which take into account the effect of changes in snow density with time on the thermal conductivity of the snow layer. By considering the dynamic nature of snow thermal conductivity, we make the snow layers consistent with the other dynamic layers of the vertical heat transfer scheme in JSBACH.  Webb et al. (2000). The time step of JSBACH is set to 30 minutes. In our simulation, we only consider the region north of 50 • N, since most of the world's permafrost soils are located in this region and it is thus possible to save a substantial amount of computation time.
We use a step-wise approach to run JSBACH into a steady-state, before we continue the simulation with a transient run leading to current climatic conditions. First, a 50-year hydrothermal spin-up is performed, where freezing and thawing of 5 water is switched off, so that only liquid soil water occurs and develops into a steady-state with climate. This is done to prevent water in the deep soil layers from freezing before the water content had time to deviate from the initial conditions, since the initial values may not necessarily be realistic. Since we do not want to include the influences of a climate change in the spin-up, these first 50 years are run with WATCH climate data randomly selected from the time period 1901 to 1930. Atmospheric CO 2 is set to a preindustrial value of 285 ppmv. Subsequently, the spin-up run is continued with a 100-year simulation, which is used 10 to generate a steady-state of soil ice, water and temperature. Therefore, freezing/thawing is switched on this time. Again, this run uses random climate data from 1901 to 1930 and preindustrial CO 2 . The hydrothermal spin-up is followed by a 5000-year simulation of the CBALANCE model, which is a simplified version of JSBACH that simulates only the slow carbon pools of soil and vegetation as a function of the NPP of the vegetation. By using the NPP output of the last 30 years of the hydrothermal spin-up and repeating it for 5000 years in CBALANCE, soil carbon is run into a steady-state. After that, the simulation is version of JSBACH and c) a configuration, where no bryophyte and lichen layer is simulated. These 3 simulations will be referred to as "Dynamic", "Constant" and "Without" througout the text. By comparing the 3 simulations we can assess the impact of the dynamic coverage and thermal properties of the bryophyte and lichen layer on soil temperature. 25 To assess how uncertainty in the parameterisation of the bryophyte and lichen PFT affects our estimated difference in soil temperature between the "Dynamic" and "Without" simulations, we run a sensitivity analysis. We test several bryophyte and lichen parameters which can affect productivity and, consequently, surface coverage and the associated thermal properties of the bryophyte and lichen layer. A detailed description of the sensitivity analysis can e found in Sect. A in the appendix.
To evaluate our modelling approach, we compare simulated bryophyte and lichen surface coverage and NPP averaged 30 over the study region to field measurements. Since large-scale observations are not available to our knowledge, we estimate "characteristic" values of surface coverage and NPP for the study region based on small-scale measurements (see Sect. B in the appendix for details).

Effects of bryophytes and lichens
Large-scale patterns of net primary productivity (NPP) and surface coverage of the bryophyte and lichen ground cover simu- The spatial pattern of NPP by bryophytes and lichens in JSBACH (Fig. 4 a) ) can be structured into areas with high productivity of over 100 g m −2 a −1 of carbon, areas of intermediate productivity of 40 to 100 g m −2 a −1 and areas of low productivity of less than 40 g m −2 a −1 . It should be noted, however, that also in areas of low productivity NPP usually exceeds 20 g m −2 a −1 .
Highly productive areas are found in North Scandinavia, the north-european part of Russia, North-West Siberia, Kamchatka, Central Canada. This spatial pattern may be explained by a combination of rainfall and available photosynthetically active radiation (PAR) on the ground (see Fig. 9 in the appendix): Regions such as North Scandinavia, North-West Russia, Kamchatka, Alaska and East Canada exhibit relatively high precipitation, but simulated above-ground vegetation is sparse enough to allow 15 for moderate levels of light on the ground. This promotes high productivity of bryophytes and lichens. South Scandinavia, parts of east-european Russia or the region south of the Hudson Bay, however, show relatively high precipitation together with little availability of PAR on the ground, which limits productivity. Large regions in North Canada and East Siberia, in turn, exhibit a high PAR availability on the ground due to sparse vegetation, but precipitation is very low there, which again limits bryophyte and lichen growth.
The bryophyte and lichen layer is excluded from JSBACH-tiles containing cropland (Sect. 2.1). Consequently, productivity of bryophytes and lichens is strongly limited by available area in regions with extensive agriculture, such as South Russia and South Canada. This explains why NPP is low in these regions in spite of favorable climatic conditions for growth.

5
The spatial distribution of simulated bryophyte and lichen surface coverage (Fig. 4 b) ) correlates in general with the spatial pattern of NPP. This is due to the uniform disturbance interval of 100 years prescribed in our simulation for the whole study region. With constant disturbance, variation in the growth of the bryophyte and lichen cover results only from variation in the organisms' NPP (see Eq. 1). However, some areas in Russia and Canada show a high surface coverage in spite of only intermediate NPP, whereas other areas show less coverage for even higher values of NPP. This is due to surface coverage being 10 plotted on a grid scale basis. If a significant part of the grid cell is covered by lakes or glaciers, coverage cannot reach 100 % although NPP might be very high. In Fig. 5 b), the annual amplitude of topsoil temperature is shown, averaged over the last 15 years of the simulations. The annual decreasing effect of the bryophyte and lichen layer on soil temperature (Fig. 5 a) ) results from a strong decrease in summer that overrules the slight increase in soil temperature in winter. Both the dynamic as well as the constant bryophyte and lichen layer show this pattern. The organic layer with constant properties, however, should actually have no overall effect on 25 annual soil temperature, since it dampens the heat flux between atmosphere and soil uniformly througout the year. Hence, to explain this outcome, the additional effect of the snow layer on vertical heat transfer in JSBACH needs to be considered: In winter, the snow layer limits heat transfer already to such an extent, that the further dampening effect of the organic layer is secondary. Interestingly, the dampening effect of the dynamic bryophyte and lichen layer is not uniform througout the year. In summer, the dynamic layer is a more efficient insulator than the constant layer, which can be explained by the additional reduc-30 tion in thermal conductivity due to low moisture content. In winter, however, the dynamic layer is less efficient in insulating the soil than the constant one, due to the occurence of ice inside the layer, which strongly increases its thermal conductivity (see also Fig. 2).
Figures 5 c) and d) display the vertical soil temperature profiles for July and January for the 7 soil layers simulated in JSBACH. In summer, the insulating effect of the bryophyte and lichen layer leads to a strong decrease in soil temperature in

Figures 6 c) and d)
show the large-scale pattern of differences in topsoil temperature and active layer thickness between the "Dynamic" and the "Without" simulation. The general pattern of soil temperature difference is similar to surface coverage, which makes sense since it is the bryophyte and lichen ground cover which reduces soil temperature. However, in mountainous 5 regions, the reducing effect of the moss and lichen cover on soil temperature seems to be stronger than in flat terrain, indicating a non-linear response of the insulating effect to the climate forcing or to soil properties. The spatial differences in soil temperature span a range of 0 to 5.7 • C.
The reduction in active layer thickness resulting from the insulating effect of the bryophyte and lichen layer is substantial, particularly at the center of the permafrost area, in Siberia and North Canada. However, also in the southern regions the reduction in soil temperature is large enough to decrease active layer thickness below the 3 m-threshold for permafrost soil. For this reason, the areal extent of permafrost soil in Fig. 6 b) is considerably larger than the coloured area in Fig. 6 d). Estimated permafrost area as defined by an active layer thickness of less than 3 m increases due to the insulating effect of the bryophyte and lichen layer by 32 % from 14.8 million km 2 to 19.5 million km 2 . This new result is much more comparable to the 22 million km 2 permafrost area reported for the Northern Hemisphere (French, 2007, p. 95). On average, active layer thickness is 5 reduced by 107 cm for the study region with a spatial range of 0 to 224 cm.
Average values over the study region are listed in Tab. 1. They show a significant impact of bryophytes and lichens on soil temperature and active layer thickness with a large spatial variation.  In this study we quantified the reducing effect of the bryophyte and lichen ground cover on soil temperature at high latitudes at the large scale. For this purpose, we estimated dynamic surface coverage and thermal properties of bryophytes and lichens with a process-based model and we integrated this model into the global land surface model JSBACH.
We estimated an average decrease in temperature of the uppermost soil layer of 2.7 K with a spatial range of 0 to 5.7 K. This is a substantial effect, it has a similar size as the projected increase in global near-surface air temperature under the RCP 4.5 warming scenario of the Intergovernmental Panel on Climate Change (Pachauri et al., 2014). Our finding is consistent with various field experiments as well as modelling studies at site level, which confirm the important role of bryophytes and lichens for reducing heat exchange between atmosphere and soil at high latitudes (Beringer et al., 2001;Gornall et al., 2007;Jorgenson et al., 2010). Our results suggest that the insulating effect of the bryophyte and lichen ground cover should be taken 10 into account in large-scale modelling studies which focus on feedbacks between permafrost soil and atmospheric CO 2 under climatic change.
Moreover, we showed that representing the dynamics of both surface coverage and thermal properties of bryophytes and lichens is crucial for estimating their insulating effect. Using a simple organic layer with constant coverage and neglecting the influence of water or ice content on thermal conductivity and heat capacity of the layer likely results in underestimating the 15 decrease in soil temperature.
Although the focus of our study is on soil temperature, the process-based bryophyte and lichen scheme in JSBACH also provides an estimate of the organisms' net primary productivity (NPP). The average simulated NPP of 49 g m −2 a −1 of carbon for the biomes "boreal forest" and "tundra" corresponds to roughly 10 % of average boreal forest NPP (Gower et al., 2001). This is a lower fraction than stated in the study of Turetsky et al. (2010), which estimate a contribution of 14 % to 58 % of 20 moss to total ecosystem NPP for various boreal forest sites. The study by Goulden and Crill (1997), however, estimates a lower contribution of around 10 % by mosses to black spruce forest NPP. The variation in the fraction of forest NPP attributed to bryophytes and lichens can be explained by differences in hydrological conditions and vegetation structure between the sites.
High values of productivity are mainly found on wetland sites, where high water and light availability sustain a productive Sphagnum-cover. JSBACH, however, is mainly designed to simulate upland soils of the boreal forest; wetlands are not yet 25 included and, consequently, the high productivity of mosses in these areas is not reflected in our estimate. Implementing a scheme for wetland hydrology and the associated additional supply of water for mosses from below would be a useful extension of JSBACH in the future.
The size of the insulating effect depended on the relations between the thermal properties of the bryophyte and lichen cover and its water or ice content. These relations are well established by field measurements and theory (see Fig. 2). Furthermore, 30 our estimate of soil temperature was sensitive to the surface coverage of bryophytes and lichens, which, in our model, largely depends on the simulated NPP of the organisms. To assess how well JSBACH is able to represent current bryophyte and lichen NPP and surface coverage at high latitudes, we compared our estimates to field measurements. However, due to large variation in the field studies and a lack in up-scaled estimates, it was difficult to define characteristic values of NPP and surface coverage for the study region (see Tab. 3 and 4). Given the considerable uncertainty in the measurements, our estimates of bryophyte and lichen NPP and surface coverage agreed well with field observations. The enhancement of JSBACH due to representation of near-surface vegetation state and functions will improve the reliability of future projections of northern ecosystem responses to environmental change as well as climate-carbon cycle feedbacks in future studies.
We compared modelled subsoil temperature and active layer thickness to observations for the region of Yakutia. Our esti-5 mated active layer thickness matches reasonably well to observations and it is improved over the previous version of JSBACH (Ekici et al., 2014). This suggests that the dynamic bryophyte and lichen layer leads to a more realistic representation of vertical heat transfer in the model. JSBACH still underestimates subsoil temperature compared to the previous model version (Ekici et al., 2014). The cold bias is most pronounced for the East-Siberian mountains where bryophyte and lichen cover is relatively low (Fig. 4 b) ). Hence, it seems likely that the reason for the cold bias is not directly related to the bryophyte and lichen layer. bryophyte and lichen cover and the environment that may affect heat exchange at the surface, which we did not consider in this study for simplicity. It has been shown by Bernier et al. (2011), for instance, that conversion of black spruce forests into lichen woodlands in Canada results in atmospheric cooling due to the higher albedo of lichens compared to forest. In JSBACH, we did not include the effect of bryophytes and lichens on surface albedo, since the organisms vary considerably in their colour.
Instead of assigning an arbitrary value for albedo to the bryophyte and lichen PFT in JSBACH, we thus assume that the albedo 20 of the PFT is similar to the soil albedo. The evaporative cooling effect of the bryophyte and lichen cover on surface temperature may be modulated by water uptake from Larch trees, which are able to root into the near-surface vegetation layer.
The bryophyte and lichen layer is represented by one single PFT in JSBACH, for reasons of consistency with vascular vegetation in JSBACH and also for computational efficiency. However, this lack of diversity may have consequences for the estimated effect of the bryophyte and lichen layer on soil temperature. Given similar climatic conditions, bryophyte and lichen 25 species may differ in their degree of water saturation and, consequently, in their thermal properties. If large regions differ in their dominant bryophyte or lichen species, this may affect our estimated patterns of soil temperature and active layer thickness. Additionally, bryophyte and lichen species differ in their thickness, while thickness within a species is relatively constant (Soudzilovskaia et al., 2013). Since the bryophyte and lichen PFT in JSBACH has a constant thickness, we may underestimate further spatial effects on the soil thermal regime.

30
Another simplifying assumption in JSBACH is the uniform disturbance interval of 100 years for the whole study region.
This value represents the average fire return interval in the boreal forest, where fire is the dominant process for disturbance (Mouillot and Field, 2005). The occurence of fire, however, may vary between regions. Fires may be much less frequent in parts of North Canada and North-East Siberia, for instance (Bonan and Shugart, 1989), which could lead to larger spatial differences in simulated bryophyte and lichen surface coverage. Under climatic change, shifts in vegetation distribution and 35 increased temperatures may lead to changes in the fire interval at high latitudes. Consequently, including a more dynamic representation of disturbance in our approach would be particularly beneficial for modelling future bryophyte and lichen cover under scenarios of climate change.

Conclusions
Here we present an new version of the global land surface model JSBACH that estimates NPP, surface coverage and dynamic 5 thermal properties of bryophytes and lichens through a process-based scheme. We apply JSBACH to quantify the impact of the bryophyte and lichen ground cover on the soil thermal regime at high latitudes. Thereby, we estimate a considerable average cooling effect of the bryophyte and lichen cover of 2.7 K (0 to 5.7 K) for the uppermost soil layer. Furthermore, we find that the strength of the cooling effect largely depends on an accurate representation of dynamic coverage and thermal properties of bryophytes and lichens. 10 These results suggest that the reducing effect of the bryophyte and lichen ground cover on soil temperature should be accounted for in studies which aim at quantifying feedbacks between permafrost soil temperature and atmospheric CO 2 due to climate change. Since our process-based approach also allows for predicting future bryophyte and lichen surface coverage at high-latitude ecosystems, a potential next step is to simulate the future impact of bryophytes and lichens on active layer thickness and permafrost extent under a transient scenario of climate change.

Appendix A: Sensitivity analysis
We perform a sensitivity analysis to assess how uncertainty in the parameter values chosen for the bryophyte and lichen PFT affects our estimates. In the following we describe several parameters which are characteristic for the bryophyte and lichen PFT and which may affect NPP, surface coverage or thermal properties of the PFT: -The parameter porosity, , of the bryophyte and lichen surface layer is set to 80 % in the model. This value is varied 5 by around 20 % in each direction, from a minimum value of 65 % to a maximum of 95 %. Porosity shows large natural variation in bryophytes and lichens, from 10 % in some lichen species (Valladares et al., 1993) to 99 % in some moss species (O'Donnell et al., 2009). We are interested, however, in uncertainty concerning the average porosity of the bryophyte and lichen ground cover for large regions. Hence, we do not vary for the full range since it is unlikely that the average porosity is close to an extreme value.

10
-The thickness z of the simulated bryophyte and lichen layer is set to 4.5 cm in the model. It is defined as the undecomposed living and dead parts of a bryophyte or lichen mat. Observed values of thallus thickness show a large variation between species, ranging from less than a millimeter to several tens of centimeters (Nash III, 1996;Bell and Hemsley, 2011). Here, we vary z from 1 to 10 cm, which is a slightly larger range than determined by Soudzilovskaia et al. (2013) for 18 bryophyte species. We do not test extreme values of z, since we are interested in bryophyte and lichen species 15 that are actually able to form a macroscopic ground cover. We do not, for instance, consider flat lichen crusts on boulder surfaces, since they are not likely to play a significant role for large-scale heat exchange between soil and atmosphere.
-The diffusivity of the water-saturated bryophyte or lichen thallus for CO 2 , D CO2,sat , affects CO 2 -uptake and, consequently, NPP and surface coverage of bryophytes and lichens in the model. We set D CO2,sat to a value of 0.01 [mol m −2 s −1 ] based on (Williams and Flanagan, 1998) and vary this value by multiplying it by the factors 0.5 and 20 2.0. We choose this form of variation since D CO2,sat shows relatively large natural variation from around 5E-4 to 2E-2 [mol m −2 s −1 ] (Williams and Flanagan, 1998;Cowan et al., 1992) and, consequently, a linear variation would not be adequate (Porada et al., 2013). It should be noted that variation in D CO2,sat represents an extension to the original bryophyte and lichen model, described in Porada et al. (2013).
-The model parameter critical water saturation, Θ crit , determines, which value of saturation is necessary for the bryophyte 25 and lichen layer to reach full metabolic activity. The increase from zero saturation and no activity to Θ crit is assumed to be linear (Porada et al., 2013). Here, Θ crit is set to 30 %, which corresponds to the lower end of the range of possible values (Porada et al., 2013) and reflects the relatively fast activation of common boreal forest floor mosses (Williams and Flanagan, 1998). Hence, we vary Θ crit by setting it closer to the upper bound, to a value of 60 %.
-The specific maintenance respiration rate R main varies over 3 orders of magnitude between different bryophyte and 30 lichen species (Porada et al., 2013). For LiBry in JSBACH, we chose an intermediate value of 1.5E-6 [mol m −2 s −1 ] and due to the large range of possible values we vary R main by multiplication by the factors 0.5 and 2.0. In the model, R main is related to photosynthetic capacity and biomass turnover rate of the bryophyte and lichen layer via the parameter "Ratio of Rubisco content to maintenance respiration", Φ RR , which describes the tradeoff between photosynthetic capacity and respiration (Porada et al., 2013). LiBry has been shown to be sensitive to this parameter, so we vary the standard value of 5 [s] by 20 % in each direction.
-In addition to photosynthetic capacity, LiBry contains 4 parameters to calculate photosynthesis, which vary between 5 different species (Porada et al., 2013): The enzyme activation energies E a,Kc and E a,Ko , which control the tempera- -Furthermore, we test how sensitive the effect of the bryophyte and lichen layer on soil temperature reacts to the disturbance interval τ D set in the model. We therefore vary the standard value of 100 years for τ D by multiplying it with the factors 0.5 and 2.0.

20
-Finally, we vary the expansion efficiency η e , which is set to 0.85, by 20 % in each direction.
For each varied parameter we run a step-wise simulation consisting of a steady-state spin-up and a transient run, as described above. For reasons of computational speed we run the model only on a single grid cell. The location of the grid cell is 55 • 30' N, 98 • 30' W, which roughly corresponds to the northern study area of the BOREAS project (e.g. Gower et al. (1997)). We then compare the difference in soil temperature of the uppermost soil level between the standard simulation and each of the 25 simulations with varied parameters to quantify the impact of parameter uncertainty on our overall estimates. Moreover, we run additional simulations with more than one varied parameter to assess their combined effect on our estimates.
Many of the tested parameters have no strong effect on NPP and surface coverage of the simulated bryophyte and lichen layer. Thus, they do not affect the difference in topsoil temperature resulting from the dynamic layer. The following parameters, however, lead to changes in the temperature difference:

30
-The porosity of the bryophyte and lichen layer affects the temperature difference in both directions. Interestingly, a reduction in porosity leads to a reduced coverage in spite of slightly increased NPP. The reason for this is that lower porosity results in a smaller specific area of the bryophyte and lichen PFT. Since net growth of cover (Eq. 1) is calculated by multiplying NPP with specific area, the effect of reduced specific area may overrule an increase in NPP. This leads to less cover expansion. The inverse effect occurs for higher porosity.
-The thickness z affects the temperature difference, to a moderate extent. The slight increase in coverage for a lower value of z mainly results from increased specific area. Comparing the parameters in Tab. 2, the slightly higher coverage 5 cannot explain by itself the increase in temperature difference. It is probably the decrease in heat capacity associated with a thinner bryophyte and lichen layer which leads to a warmer and, consequently, drier layer in summer which acts as a more efficient insulator. The opposite then happens for a thicker layer. This means that the reducing effect of low moisture content on thermal conductivity of the bryophyte and lichen layer overrules the increasing effect of higher thickness in the model.

10
-Reducing the expansion efficiency η e also reduces surface coverage and, consequently, the soil temperature difference. However, the model is not very sensitive to this parameter, η e would have to be quite low to significantly affect our estimates.
-Increasing critical water saturation Θ crit has a relatively strong decreasing effect on the temperature difference, since the associated slower activation of bryophyte and lichen photosynthesis significantly reduces NPP and surface coverage.

15
-The specific maintenance respiration rate R main has the strongest influence on the temperature difference compared to the other parameters, but only in one direction. A doubling of R main significantly reduces the temperature difference, whereas halving R main does not have any effect. The reason for this is that our standard value of R main is close to optimal, meaning that NPP decreases in both directions. The associated biomass turnover, however, increases linearly with R main and it has a strong reducing impact on net growth and cover expansion. This means that low turnover 20 compensates for low NPP at low R main , but high turnover has an additional reducing impact at high R main . The "Ratio of Rubisco content to maintenance respiration", Φ RR , has only a minor impact on the temperature difference.
-Halving the molar carboxylation rate of Rubisco, V M,C , reduces significantly the temperature difference due to strongly decreased NPP and coverage. Doubling V M,C , however, does not have any effect on temperature due to the nonlinear response of coverage on increased NPP. The activation energy E a,Kc has a moderate reducing impact on NPP, coverage 25 and thus on the temperature difference.
-Reducing the optimum temperature of photosynthesis, T opt , by 5 • C results in significantly lower NPP, coverage and temperature difference, while increasing T opt has no significant effect. The reason for the asymmetric response of NPP to T opt is the exponential dependence of respiration on the, in this case, mostly negative difference between surface temperature and T opt . Since T opt is already relatively high, a further increase leads only to a small reduction in respiration. To assess the effect of combined varied parameters, we select four parameters which significantly affect the estimated difference in topsoil temperature between the "Dynamic" and "Without" simulations (see Tab. 2). The outcome of this analysis is shown in Fig. 7. It can be seen that the combined effect of parameters is not additive, but becomes weaker with each additional parameter.

Appendix B: Model evaluation
To evaluate our modelling approach, we compare the JSBACH estimates of bryophyte and lichen surface coverage and NPP to field observations. However, to our knowledge, there are no field measurements of bryophyte and lichen surface coverage 5 and NPP which cover the global scale. Therefore, we compile a list of available small-scale observations from various studies, which are largely taken from overviews in Turetsky et al. (2010) and Bona et al. (2013). We do not attempt to create a comprehensive review of measurements of surface coverage and NPP in high-latitude regions. However, our list is sufficient to establish "characteristic" values of surface coverage and NPP for the region north of 50 • N. Thereby, we constrain our analysis to the biomes "boreal forest" and "tundra" (after Olson et al. (2001), see also Fig. 9 c) ). Since JSBACH is not primariy measurements both on peatland and upland sites. In this case, we do not consider values from peatlands. In many cases, the original studies provide NPP measurements in grams of biomass. We convert these values to grams of carbon using a factor of 0.45 (Bauer et al., 2009).
We compare the characteristic values of measured NPP to the average bryophyte and lichen surface coverage and NPP simulated by JSBACH, also constrained to boreal forest and tundra. Furthermore, we compare simulated soil temperature and 5 active layer thickness to large-scale, observation-based maps for the region of Yakutia. Table 3 shows studies which provide measurements of bryophyte and lichen NPP for the biomes "boreal forest" and "tundra".
To obtain a characteristic value for NPP, we calculate the median NPP of each study and from this set of values we take again the median, which is then shown in the last row of Tab. 3 together with the range of the median values from all studies. Given the large variation in these median values of two orders of magnitude, it is difficult to constrain NPP in the study region based on field measurements. The average value probably lies somewhere between 10 and 100 g m −2 a −1 of carbon.
In Tab. 4 several studies are listed which measure surface coverage of bryophytes and lichens for boreal forest and tundra.
Same as with NPP, we show in the last row the range and median of all studies' median values to obtain a characteristic surface coverage for the region. The range of observed surface coverage is large, but half of the (rounded) values lie between 0.6 and 0.8. We want to point out that the study by Rapalee et al. (2001) provides a large-scale estimate of surface coverage, based on 15 remote sensing. They cover the whole BOREAS region which has an area of approximately 500000 km 2 and arrive at a moss surface coverage of 0.57. This value agrees well with our estimate, given that the BOREAS region has a slightly lower than average simulated bryophyte and lichen coverage (Fig. 4). Figure 8 compares the spatial patterns of subsoil temperature and active layer thickness estimated by JSBACH to observations from Beer et al. (2013) for the region of Yakutia. JSBACH underestimates subsoil temperature to a similar extent as the 20 version presented in Ekici et al. (2014). Potential reasons for the mismatch are discussed in Ekici et al. (2014) and include a reduced spatial heterogeneity of the observational map or biases in the representation of climate, soil properties or snow depth in JSBACH. The estimate of active layer thickness, however, is improved in general compared to the previous version of JSBACH, with the exception of the East-Siberian mountain range. The underestimation of active layer thickness there likely follows from the strong underestimation of soil temperature in this area. The overestimation of active layer thickness in the Table 3. Studies which measure NPP in the biomes "boreal forest" and "tundra". Median and range of all values listed in a study are shown, "#" stands for the number of values. NPP is in g m −2 a −1 of carbon. "BOREAS N or S" stands for the northern and southern site of the BOREAS project (e.g. Gower et al. (1997)). BOREAS N is located in Manitoba, Canada and BOREAS S is located in Saskatchewan, Canada.
The last row shows median and range of the values in the first column.

Median
Range  (Skre and Oechel, 1979) Central Alaska Table 4. Studies which measure surface coverage in the biomes "boreal forest" and "tundra". Median and range of all values listed in a study are shown, "#" stands for the number of values. The BOREAS N study area is located in Manitoba, Canada and BOREAS S is located in Saskatchewan, Canada (Gower et al., 1997). The last row shows median and range of the values in the first column. Spatial patterns of a) precipitation and b) photosynthetically active radiation (PAR) on the ground. Note that PAR on ground is a variable of LiBry in JSBACH, it is therefore only larger than zero where bryophytes and lichens exist in the model. c) shows the biome mask based on Olson et al. (2001) which is used to constrain our estimates to "boreal forest" and "tundra".

Median
Author contributions. P. Porada developed the model code and prepared the manuscript with contributions from all co-authors, P. Porada and C. Beer designed the simulation setup and performed the simulations.