Modelling the diurnal and seasonal dynamics of soil CO 2 exchange in a semiarid ecosystem with 1 high plant-interspace heterogeneity 2

12 This study represents a first attempt to model the diurnal and seasonal dynamics of soil CO2 exchange 13 (FS) in a dryland ecosystem with a high plant-interspace heterogeneity. The modelling used an 14 integrated process-based approach, in which the CO2 production, transport and surface exchanges (e.g. 15 biocrust photosynthesis, respiration and photodegradation) are considered simultaneously. The model 16 was parameterized and validated with multivariate data measured during year 2013-2014 in a semiarid 17 shrubland ecosystem in Yanchi, northwestern China. We also investigated the sensitivity of simulated 18 FS to a set of stand-specific parameters and investigated the relative contribution of different flux 19 components. The model explained reasonably well the two-year dynamics of FS measured from a non20 crusted and two lichen-crusted plots. Simulations showed that the temporal pattern of FS could deviate 21 from that of the total CO2 production from rooting-zone soil. Such deviations could be explained by 22 the variations of CO2 dissolution and the CO2 exchanges of biocrust during wetting-drying cycles, and 23 the root uptake and transport of dissolved CO2. Moreover, the FS was spatially sensitive to the plant24 interspace differences and the variations in root biomass, soil organic matter and pH. These results 25 emphasized that, the processes beyond autotrophic and heterotrophic respirations and the 26 heterogeneities of soil at plant-interspace can strongly affect the FS dynamics and their climatic 27 sensitivities. Such variability should be carefully considered in extrapolation of findings from 28 chamber to ecosystem level and from seasonal to inter-annual scales. Based on this work, our model 29 can serve as a useful tool to simulate FS dynamics in dryland ecosystems. 30


Introduction
CO 2 exchange between soil and atmosphere constitutes a major carbon (C) loss from terrestrial ecosystems (Raich et al., 2002;Giardina et al., 2014).It also plays an important role in the feedbacks between the global carbon cycle and climate change (Rustad et al., 2000;Giardina et al., 2014;Karhu et al., 2014).Arid and semiarid (dryland) ecosystems cover over 40 % of land surface and contribute notably to interannual variations of terrestrial C sink (Poulter et al., 2014).However, the contribution of soil CO 2 flux (F S ) from those ecosystems to the global C budget is less studied (Castillo-Monroy et al., 2011;Gao et al., 2012;Jia et al., 2014).The temperature dependency of biological CO 2 production (i.e.autotrophic respiration and heterotrophic respiration) serves a conventional basis for F S modelling in many terrestrial ecosystems (Raich and Tufekciogul, 2000;Ryan and Law, 2005;Song et al., 2015).Soil CO 2 flux of dryland ecosystems is also widely interpreted using temperature-response functions modified by other environmental constraints, e.g.soil water content, abundance of substrates and microbial activities (Curiel Yuste et al., 2007; W. Wang et al., 2014;B. Wang et al., 2014B. Wang et al., , 2015)).
Although many empirical studies have explained the dynamics of soil CO 2 flux in specified space-time, their lack of mechanistic descriptions represents a major difficulty in extrapolation under changing environmental conditions (Fan et al., 2015).Soil CO 2 flux is a "bulk" exchange that comprises two main sets of processes, i.e. the CO 2 production and transport (Fang and Moncrieff, 1999;Fan et al., 2015).Hence, models considering only autotrophic and heterotrophic respiration often fail to account for the observed F S dynamics (Austin and Vivanco, 2006).Gas-transport processes are im-Published by Copernicus Publications on behalf of the European Geosciences Union.
portant mechanisms regulating the magnitudes and hysteretic features of soil CO 2 fluxes (Ma et al., 2013).A substantial fraction of respired CO 2 may be transported to the atmosphere via xylem and may not be measured by techniques like soil reparation chambers (Bloemen et al., 2013(Bloemen et al., , 2016)).During wet periods, soil CO 2 efflux could decrease significantly by water clogging of soil pores, which restricts the diffusion of O 2 and CO 2 gases (Šimunek and Suarez, 1993;Fang and Moncrieff, 1999).In dryland soils of high salinity/alkalinity, CO 2 transport and the water cycle are tightly coupled, as large inorganic C fluxes can be driven solely by dissolution and infiltration of CO 2 and carbonates (Buysse et al., 2013;Ma et al., 2013;Fa et al., 2014).Such inorganic processes may introduce fluctuations not only to hourly or diurnal soil CO 2 effluxes (e.g.Emmerich, 2003;Xie et al., 2009;Buysse et al., 2013) but also to terrestrial CO 2 sinks at much broader spatiotemporal scales (Schlesinger, 2009;Li et al., 2015).
Key processes contributing to CO 2 production in dryland soils also extend beyond autotrophic and heterotrophic respiration.Although biocrust organisms (lichens, mosses, bacteria, fungi and microfauna) mainly inhabit in the top few centimetres of a soil profile, they constitute up to 70 % of biomes in interspace areas (Belnap et al., 2003).These communities are able to uptake C from atmosphere (Belnap et al., 2003;Castillo-Monroy et al., 2011;Maestre et al., 2013), leading to greater concentrations of organic matter in the soil crust than in the layer underneath (Ciais et al., 2013).Although crust organisms could be inactive under stresses (e.g.drought; Green and Proctor, 2016), their photosynthetic potentials could be large (Zaady et al., 2000;Lange, 2003), even comparable to temperate forests with closed canopies (e.g.Zaady et al., 2000).The C uptakes of biocrusts are highly sensitive to stresses like droughts, thermal extremes and excessive ultraviolet radiation (Pointing and Belnap, 2012).Such variations can readily alter crusted soils between considerable CO 2 sinks and sources within a few hours (e.g.Bowling et al., 2011;Feng et al., 2014).In addition, the accumulation of debris from crust and canopy fuel photodegradation represents an important abiotic C loss in arid conditions aside from biotic decompositions (e.g.Austin and Vivanco, 2006;Throop and Archer, 2009).Photodegradation is likely to dominate the mineralization during dry daytime periods, when radiation is strong and microbial activities are prohibited by low moisture content and high temperatures (e.g.Gliksman et al., 2016).On an annual basis, photodegradation could consume more than 10 % of soil organic matter (SOM) at the surface (e.g.Austin and Vivanco, 2006;Henry et al., 2008;Brandt et al., 2010), even for the substrates (e.g.lignin) that are difficult to degrade via biotic pathways (Henry et al., 2008).
The influences of the multiple C processes (i.e.autotrophic and heterotrophic respiration, C exchanges by biocrust organisms, inorganic C fluxes and photodegradation) on soil CO 2 exchanges are highly overlapped and tightly related to water-energy conditions.In dryland ecosystems, patchy veg-etation and large fractions of interspace area are common features (Domingo et al., 2000), and the water-thermal conditions can vary considerably from plant-covered areas to interspace within even a few metres (Rodríguez-Iturbe et al., 2001;Caylor et al., 2008;Ma et al., 2011).The water-energy dynamics at the different surfaces are linked by multiple advection processes both above and below the ground (Gong et al., 2016).Due to the complexity of water-energy processes, there may exist high non-linearity of water-thermal responses to climatic variabilities (e.g.Phillips et al., 2011;Barron-Gafford et al., 2013).This will also complicate the C responses and consequently affect the relationships between CO 2 fluxes and environmental controls (e.g.Jarvis et al., 2007;Song et al., 2015).
Global climate change is expected to increase annual mean air temperatures considerably and alter precipitation regimes (Donat et al., 2016).Understanding the responses of dryland ecosystems to such changes requires mechanistic models that integrate the multiple biotic and abiotic processes in soil C cycling.So far, only a few models have coupled the biotic CO 2 production with the transport of gases and heat (Šimunek and Suarez, 1993;Fang and Moncrieff, 1999;Phillips et al., 2011;Ma et al., 2013;Fan et al., 2015).Nevertheless, none of those models have described the heterogeneous water-energy conditions in the soil-plant-atmosphere continuum (SPAC) nor the unconventional C fluxes such as C uptake by biocrusts and photodegradation, despite the importance of these processes to arid and semiarid environments.Perhaps models by Porada et al. (2013) and Kinast et al. (2016) represent the few existing works in this sense.However, both models focus on the patterns at the regional scale with very simplified ecosystem processes.In this study, we aim to (i) investigate the roles of soil CO 2 components in regulating soil CO 2 effluxes in a dryland ecosystem using a process-based modelling approach and (ii) estimate the plant-interspace differences in the componential C processes.

Model overview
The process-based model was build based on a semiarid shrubland ecosystem located at the southern edge of Mu Us Desert (37 • 42 1 N, 107 • 13 7 E; 1560 m above sea level; Fig. 1a), Ningxia, China (see Wang et al., 2014Wang et al., , 2015)).The long-term  mean temperature is 8.1 • C, and the mean annual precipitation is 287 mm, most of which falls from July to September (Jia et al., 2014).The radiation and evaporation demands are high in this area.The annual incoming shortwave radiation is 1.4 × 10 5 J cm −2 and the annual potential evaporation is 2024 mm.The vegetation is dominated by scattered crowns of Artemisia ordosica (Fig. 1b).The soil is highly alkaline (pH of 8.2).Biocrust (mainly lichens and algae) covers about 40 % of interspace soil (Fig. 1c-e).The thickness of the crust layer was 0.5-2.5 cm (Gong et al., 2016).
In the modelling, the structure of ecosystem was considered as replications of representative land units (RLUs; Fig. 1f; Gong et al., 2016), which comprise the area covered by shrubs and the surrounding soil (interspace).Vertically, the model simulates the C flows across soil profile and the water-energy transport from the lower boundary of rooting zone to a reference height in the boundary atmosphere.Horizontally, the SPAC processes at the plant-covered and interspace areas are differentiated but related via advection and diffusion flows, as driven by the gradients of temperature, water potential and gas concentration.The mineralization, uptake and transport of soil C and nitrogen (N) are also regulated by water-energy conditions.
Figure 2 shows the framework of key processes and variables included in the F S modelling.The model includes a set of submodels which describe (i) CO 2 dissolution, transport and efflux; (ii) autotrophic and heterotrophic CO 2 production in the soil profile; (iii) CO 2 uptake and emission by biocrust; (iv) surface energy balance and the soil temperature profile; and (v) soil hydrology and water balance.These submodels are linked by multiple feedbacks to represent the coupling of C, water, vapour and energy transportation in the ecosystem.Submodels (iv)-(v) have been developed and described in detail in our previous work (Gong et al., 2016), which focused on (i) introducing the plant-interspace heterogeneity into water-energy modelling and (ii) investigating the influences of such heterogeneity on the ecosystem water-energy budgets for a dryland ecosystem.Gong et al. (2016) also validated the model in regard to the diurnal to seasonal dynamics of radiation balance, surface energy balance, soil temperature and moisture content in the footprint area of a eddycovariance (EC) site (for details of measurement, see Jia et al., 2014).In this work, we therefore focused on the development of submodels (i)-(iii) and the model validation using F S data measured by multiple automatic respiration chambers from crust-covered and non-crusted soils.Based on the validated model, a series of sensitivity analyses was carried out addressing (i) how the soil CO 2 components regulate F S in the studied ecosystem and (ii) how the plant-interspace heterogeneities influence the componential C processes and F S .

Submodel (i): CO 2 transport, dissolution and efflux
For soil fraction x (see Fig. 1f for RLU settings), CO 2 exchange (F S , upward positive) was the sum of CO 2 uptake by biocrust (F C t ), photodegradation (F P ) and the emission from soil under the biocrust layer (F T ): where F C t is the net balance between biocrust photosynthesis (P Ct ) and respiration (R Ct ), and F C t = P Ct − R Ct (see Sect. 2.2.3).F T was modelled based on the mass-balance functions developed by Fang and Moncrieff (1999), which combined major transport processes in both gaseous and liquid phases.We expanded the original function from one dimension to two dimensions.For the soil layer (x, i) and time step t, the CO 2 concentration and C flows were calculated as follows: where superscripts v and h denote the vertical and horizontal directions, respectively (see also in Gong et al., 2016); C is the total CO 2 content; F dg and F dw are the CO 2 flows due to diffusion/dispersion via the gaseous and liquid phases; F ag and F aw are the flows in gaseous and liquid phases due to gas convection and water movement; and S is the net CO 2 sink of the layer.The calculation schemes of F dg , F dw , F ag and F aw have been described in detail by Fang and Moncrieff (1999).F T is the total exchange of gaseous CO 2 between the surface and topmost layer: where E S x,1 is the soil evaporation at section x (see Eq. 17) in Gong et al. (2016); C w is the equivalent CO 2 concentrations in the soil solution.For layer (x, i), C w is linked to the gaseous CO 2 concentration (C g ): where V is the total porosity, and θ is soil water content.C g and C w were further related via the dissolutiondissociation balance of CO 2 in the soil solution, following Fang and Moncrieff (1999) and Ma et al. (2013): where P C is the partial pressure of CO 2 in pore air; K H is Henry's law constant; K 0 , K 1 and K 2 are the equilibrium coefficients of dissolution, the first-and the second-order dissociation reactions of carbonic acid, respectively (for details, see Fang and Moncrieff, 1999).The equilibrium [H + ] was determined by soil pH and the coefficients K H , K 0 , K 1 and K 2 , which were functions of temperature in each soil layer (Fang and Moncrieff, 1999).C w was calculated as the sum of CO aq 2 , H 2 CO 3 , HCO − 3 and CO 2− 3 .2.2.2 Submodel (ii): autotrophic and heterotrophic CO 2 production along the soil profile For soil layer (x, i), S x,i (Eq.2) was calculated as the sum of autotrophic and heterotrophic CO 2 production, and the dissolved CO 2 was removed with the water taken up by roots: where E is the transpiration uptake of water (Gong et al., 2016); R s is the CO 2 production by heterotrophic SOM decomposition; R a is the autotrophic respiration of the rhizosphere, which comprises maintenance respiration (R m ) and growth respiration (R g ): To simulate R s , we simplified the pool-type model of Gong et al. (2013Gong et al. ( , 2014)), which originated from Smith et al. (2010) for simulating coupled C and N cycling in organic soils.
The SOM pool in each soil layer was divided into debris (M deb , i.e. litter from roots and biocrust), microbes (M mic ) and humus (M hum ), which are different in biochemical recalcitrance and N content.During decay, mineralized masses transfer from M deb and M mic to a more resistant form (i.e.M hum ), leading to a decrease in labiality (e.g.Li et al., 1992).The mineralization of organic C followed first-order kinet-ics and was constrained by multiple environmental multipliers, including temperature, water content and oxygen content (Šimunek and Suarez, 1993;Fang and Moncrieff, 1999): where superscript r denotes the type of SOM pool (r = 1 for M deb , r = 2 for M mic and r = 3 for M hum ); m is mineralized SOM during time step dt; k is the decomposition constant; dt is the time step; f (T s x,i ), f (θ x,i ) and f (O x,i ) are multiplier terms regarding the temperature, water content and oxygen restrictions, respectively.f (O x,i ) was calculated following Šimunek and Suarez (1993).f (T s x,i ) and f (θ x,i ) were reparameterized with respect to the site-specific conditions of plants and soil (see Sect. 2.4.3).The CO 2 production from mineralization was further regulated by the N starvation of microbes following Smith et al. (2010): where r E is the gas production rate (r E [0, 1]), and (1−r E ) is the proportion of organic matter passed to downstream SOM pools.The evolution of each SOM pool was calculated as below: where A is the SOM input rate (A = 0 for M mic and M hum ); superscript r − 1 denotes the source SOM pools.R m x,i was calculated in a way similar to R s x,i (e.g.Chen et al., 1999;Fang and Moncrieff, 1999).R g x,i was calculated as a fraction of photosynthetic assimilates, following Chen et al. (1999): where M R is the root biomass; k R is the specific respiration rate of roots; k g is the fraction of photosynthetic assimilation consumed by growth respiration; f r x,i is the mass fraction of roots in the soil layer (x, i).P g is the photosynthesis rate of plants.P g was estimated using a modified Farquhar's leaf biochemical model (see Chen et al., 1999).This model simulates the photosynthesis based on biochemical parameters (i.e. the maximum carboxylation velocity, V max and maximum rate of electron transport, J max ), foliage temperature (T c ) and stomatal conductance (g s ).The values of V max and J max were obtained from in situ measurements from the site (Jia et al., unpublished).T c and g s were given in the energy balance submodel, which was detailed in Gong et al. (2016).N content bonded in SOM is mineralized and released to soil layers simultaneously with decay.The abundance of mineral N (i.e.NH + 4 and NO − 3 ) regulates the growth of microbial biomass and r E following Smith et al. (2010) and Gong et al. (2014).Key processes governing the dynamics of mineral N pools included nitrification-denitrification (Smith et al., 2010), solvent transport with water flows (Gong et al., 2014) and the N uptake by the root system.However, the plant growth was not modelled in this work.Instead, N upt was calculated using the steady-state model of Yanai (1994), based on the transpiration rate, surface area of fine roots and the diffusion of solvents from pore space to root surface: where r o is the fine root diameter; L is the root length and 2π r o L is the surface area of fine roots; α is the nutrient absorbing power, which denotes the saturation degree of solute uptake system (α ∈ [0, 1]); C o is the concentration of solvents at the root surface and is a function of bulk concentration of mineral N (N min ), inward radial velocity of water at root surface (v o = E (2π r o L)) and saturation absorbing power α.Further details for calculating α and C o can be found in the work of Yanai (1994).

Submodel (iii): CO 2 exchange of biocrust and photodegradation
Biocrusts are vertically layered systems that comprise the top crust (or bio-rich layer) and underlying subcrust (or biopoor layer), which are different in microstructure, microbial communities and C functioning (Garcia-Pichel et al., 2016;Raanan et al., 2016).Top crust is usually a few millimetres thick, which allows the penetration of light and the development of photosynthetic microbes (Garcia-Pichel et al., 2016).On the other hand, the subcrust has little photosynthetic activity.Here, we focused mainly on describing the C exchanges in the top crust but assumed the C processes in the subcrust were similar to those in the underneath soil.We developed the following functions to describe the C fixation and mass balance in the top crust: where P Ct is the bulk photosynthesis rate, and R Ct is the bulk respiration rate.P C and R C were further modelled as follows: where α C is the apparent quantum yield; P C m is the maximal rate of photosynthesis and was a function of moisture content (θ Ct ) and temperature (T Ct ) in the top crust; A PAR is the photosynthetically active radiation (PAR); M Ct is the total C in the SOM of top crust; k cr is the respiration coefficient; f (θ Ct ) and f (T Ct ) are water and temperature multipliers.Here, we assumed zero photosynthesis rate for subcrust.The heterotrophic respiration (R C s ) was calculated following Eq.( 11), based on the C storages (M r x,1 ), temperature and moisture content of the crust layer (i.e.T s x,1 and θ x,1 ; see Eqs. 29 and 14 in Gong et al., 2016).
To consider different C losses and exchanges, and to calculate the C balance in top crust and subcrust, respectively, we considered the following matters.R Ct includes the respiration from both autotrophic (M CA ) and heterotrophic (M CH ) pools.When autotrophic organisms die, SOMs pass from M CA to M CH and influence the turnover processes.A variety of top crust organisms can reach into subcrust (e.g. through rhizines; Aguilar et al., 2009) and export litter there.When the surface is gradually covered by deposits, top crust organisms tend to move upward and recolonize at the new surface (e.g.Garcia-Pichel and Pringault, 2001;Jia et al., 2008), leaving old materials buried into the subcrust (Felde et al., 2014).On the other hand, the debris left to soil surface is exposed to photodegradation.Based on the above information, the C balance in top crust and subcrust was calculated as follows, assuming the partitioning of respiration between autotrophic and heterotrophic pools was proportional to their fractions: where k m is the rate of C transfer (e.g.mortality) from autotrophic pool to heterotrophic pool; k b is the rate of C transfer (e.g.burying) from top crust to subcrust; F P is the loss of SOM due to photodegradation.Photodegradation tends to decrease surface litter masses in a near-linear fashion with the time of exposure (Austin and Vivanco, 2006;Vanderbilt et al., 2008).Considering the diurnal and seasonal variations of radiation, F P was calculated as a function of surface SOM mass and solar radiation: where Rad x is the incident shortwave radiation at surface x (Gong et al., 2016); M surf is the surface litter mass; and k p is the photodegradation coefficient.

Micrometeorological and soil CO 2 efflux measurements
Meteorological variables were measured every 10 s and aggregated to half-hourly resolution during 2013-2014.The factors measured included the incoming and outgoing irradiance (PAR-LITE, Kipp and Zonen, the Netherlands), PAR (PAR-LITE, Kipp and Zonen, the Netherlands), air temperature and relative humidity (HMP155A, Vaisala, Finland).Rainfall was measured with a tipping-bucket rain gauge (TE525WS, Campbell Scientific Inc., USA) mounted at a nearby site (1 km away; see B. Wang et al., 2014).The seasonal trends of the measured T a and P can be found in Jia et al. (2016).No surface runoffs were observed at the site, indicating the horizontal redistribution of rainfall was mainly through subsurface flows.
Continuous measurements of F S were conducted using an automated soil respiration system (model LI-8100A fitted with a LI-8150 multiplexer, LI-COR, Nebraska, USA).The system was on a fixed sand dune of typical size (B.Wang et al., 2014), which was located about 1.5 km south of the EC tower described in Gong et al. (2016).Three collars (20.3 cm in diameter and 10 cm in height, of which 7 cm were inserted into the soil) were installed, on average, at 3 m spacing in March 2012.One collar (C1; see Fig. 1c) was set on a bare-soil microsite with no presence of biocrust.Two other chambers (C2, see Fig. 1d; C3, see Fig. 1e) were set on lichen-crusted soils.F S was measured hourly from C1 and C2 by opaque chambers, whereas it was measured by transparent chamber from C3 to include the photosynthesis and photodegradation.Litter from the shrub canopies was cleared from the collars during weekly maintenance.Hourly T s and θ at 10 cm depth were measured outside each chamber using the 8150-203 soil temperature sensor and ECH2O soil moisture sensor (LI-COR, Nebraska, USA), respectively.Root biomass was sampled near each collar (within 0.5 m) in July 2012, using a soil corer (5 cm in diameter) to a depth of 25 cm.The samples were mixed and sieved sequentially through 1, 0.5 and 0.25 mm meshes, and the living roots were picked by hand.The comparison of the three microsites is shown in Table 1.Methods used in data processing and quality control have been described earlier in detail (see Wang et al., 2014Wang et al., , 2015)).The quality control led to gaps of 10-13 % in the F S dataset.

Parameterization of vegetation and soil texture
The parameterization schemes supporting the simulations of energy balance and soil hydrology in submodels (i)-(v) have been described previously in detail by Gong et al. (2016).As the water-energy budget is sensitive to vegetation (i.e.canopy size, density and leaf area) and soil hydraulic properties (see Gong et al., 2016), we hereby re-estimated these parameters for the F S site.Measurements based on four 5 m × 5 m plots showed that the crown diameter D (86 ± 40 cm) and height H (47 ± 20 cm) at this site were similar to those measured from the eddy-covariance (EC) footprint by Gong et al. (2016).However, the shrub density was 50 % greater, leading to higher shrub coverage (42 %), shorter spacing distance L (40.2 cm) and greater foliage area.On the other hand, the subsoil at the F S site is sandy and much coarser than that at the EC footprint.Therefore, we collected 12 soil cores from 10 cm depth and measured saturated water content (θ sat ), bulk density and residual water content (θ r ) from each sample.Then, the samples were saturated, and covered and drained by gravity.We measured the water content after Table 2. Parameters for soil water retention and C turnover.

Parameter Equation Unit
Value 2 and 24 h draining, which roughly represented the matrix capillary water content (10 kPa) and field capacity (33 kPa) (Armer, 2011).The shape parameters n and α h (see Eq. 26 in Gong et al., 2016) for the water-retention function were estimated from these values (Table 2).

Parameterization of soil C and N pools
The sizes and quality of soil C pools were parameterized based on a set of previous studies.The total soil organic carbon (SOC) in the root-zone soil (i.e. 60 cm depth, bulk density of 1.6 g cm −3 ) was set to 1200 g m −2 , based on the values reported from previous studies in Yanchi area (e.g.Qi et al., 2002;Chen and Duan, 2009;Zhang and Hou, 2012;Liu et al., 2015;Lai et al., 2016).The mass fraction of the resistant SOM pool (M hum ) was set to 40-50 % of the total SOM, following work of Lai et al. (2016).The vertical distribution of the SOM pools was described following Shi et al. (2013).At the ecosystem level, the aboveground biomass was linearly related to the crown projection area (M S = 0.2917 × π (0.5D) 2 ; see Zhang et al., 2008).The total root biomass was then calculated as proportional to the aboveground biomass using a root-shoot ratio of 0.47 (M R = 0.47M S ; Xiao et al., 2005).The vertical profile of root biomass was parameterized as decreasing exponentially with depth, using the depth profile reported by Lai et al. (2016).In the horizontal direction, root biomass was set to decrease linearly with the distance from the centre of a shrub crown (Zhang et al., 2008).The N content was parameterized following the measurement of Wang et al. (2015).
Based on the above settings, the specific decomposition rate of debris was estimated from the litterbag experiment done by Lai et al. (2016), which showed a 16 % decrease in the mass of fine-root litter during a 7-month period in 2013 at the Yanchi site.The photodegradation coefficient (k p ) was set to 0.23 yr −1 , which was the mass-loss rate reported by Austin and Vivanco (2006).M surf was set to 33 % of M CH in top crust, assuming the depth of light penetration was about 2 mm and C concentration was homogeneous in top crust.The surface litter from canopy was not considered in this modelling, as the plant litter was cleaned from the collars during weekly maintenance.The specific respiration rate of roots (k R ), however, could be much greater during vegetative growing stage than other periods, e.g. at the defoliation stage (Fu et al., 2002;Wang et al., 2015).Here, we linked k R to the development of foliage in modelling using the approach of Curiel Yuste et al. (2004): where k R0 is the "base" respiration rate (Table 2); L l is the green leaf area, which is a function of the Julian day (Gong et al., 2016); L max is the maximum L l ; n R is the maximum percentage of variability and is set to 100 %.

Parameterization of soil CO 2 production
Based on the empirical study of B. Wang et al. (2014), the steady-state sensitivity of CO 2 production to soil temperature and water content (i.e.f (T s ) f (θ ); Eq. 11) can be described as a logistic-power function: where a, b and c are empirical parameters.This function represents the long-term water-thermal sensitivity of CO 2 production over the growing seasonal, yielding an apparent temperature sensitivity Q 10 of 1.5 for the emitted CO 2 (B.Wang et al., 2014).However, this could underestimate the shortterm sensitivities of CO 2 production.The apparent Q 10 could be much greater at the diurnal level than at the seasonal level (B.Wang et al., 2014).In this work, we firstly calculated the base sensitivity using the long-term scheme (Eq.26) with a 1-day moving average of water-thermal conditions.Then, the deviation of hourly sensitivity from the base condition was adjusted by the short-term Q 10 : Q 10 (θ short ) = 18 010θ 3.721 short + 1.604, where T s short and θ short are the 1-day moving averages of T s and θ , respectively; Q 10 (T s ) and Q 10 (θ ) are the adjustment functions for short-term apparent Q 10 , regarding the shortterm T s and θ .Further non-linearity of soil respiration responses refers to the rain-pulse effect (or the "Birch effect"; Jarvis et al., 2007) that respiration pulses triggered by rewetting can be orders of magnitude greater than the value before the rain event (Xu et al., 2004;Sponseller, 2007).Such a response could be very rapid (e.g.within 1 h to 1 day, Rey et al., 2005) and sensitive to even minor rainfall.It also seems that the size and duration of a respiration pulse not only depend on the precipitation size but also on the moisture conditions prior to the rainfall (Xu et al., 2004;Rey et al., 2005;Evans and Wallenstein, 2011).As numerical descriptions of such an effect remain unavailable at the moment, we simply multiplied Eq. ( 26) to a rain-pulse coefficient (f pulse ): where θ 72 h is the 3-day moving average of soil moisture content; n p is a shape parameter and set to 2 in this study.θ 72 h is the 72 h moving average of θ.For tests on model sensitivities to different parameterizations of f pulse , see Sect.2.5.3.

Parameterization of biocrust photosynthesis and respiration
In submodel (iii), Eqs. ( 17)-( 19) were parameterized based on the experiment of Feng et al. (2014).In the experiment, 50 lichen (top crust) samples of 0.5-0.7 cm thickness (100 % coverage; average C content of 1048 µmol C cm −3 ) were collected from a 20 m × 20 m area.The samples were wetted and incubated under controlled T Ct (i.e.35, 27, 20, 15 and 10 • C).These samples were divided into two groups to measure the net primary productivity (NPP) and dark respiration (R d ) separately.Gas exchanges and light response curve for each crust sample were measured using an LI-6400 infrared gas analyser equipped with an LI-6400-17 chamber and an LI-6400-18 light source (LI-COR, Nebraska, USA).Measurements were taken at ambient CO 2 values of 385 ± 35 ppm.Saturated top crust samples were placed in a round tray and moved to the chamber.CO 2 exchange was measured during the drying of samples until the CO 2 flux diminished.During drying, θ Ct was measured every 20 min.
For more details, see Feng et al. (2014).
Fitting measured R d to T Ct and θ Ct (see Fig. 3a) was based on the Matlab ® (2012a) curve-fitting tool.The obtained multipliers in Eq. ( 19) are as follows: where Q Ct , a R C , b R C and c R C are the fitted shape parameters (Table 2).The parameterized Eq. ( 19) was then used to simulate the R d for the NPP samples, based on the correspondent T Ct and θ Ct from each measurement.P C m was determined by subtracting the simulated respiration rate from the NPP measured under light-saturated conditions.Then, P C m was fitted to T Ct and θ Ct in the Matlab ® (2012a) curve-fitting tool using the following equations (Fig. 3b): where a P t , b P t , c P t , d P t , a P w , b P w , c P w and d P w are fitted shape parameters (Table 2).
It should be addressed that T Ct and θ Ct could change more rapidly than the mean conditions of the crust (i.e.T s x,1 and θ x,1 ).In this work, T Ct was calculated from the surface temperature (T x ; see Eq. 13 in Gong et al., 2016) and T s x,1 by linear interpolation.The calculation of θ Ct , on the other hand, depended on the drying-rewetting cycle.During drying phases, θ Ct was interpolated linearly from θ x,1 and surface moisture content (θ x ), whereas during wetting phases the mass balance of water input P and evaporation loss (E s see Eq. 17 in Gong et al., 2016) was considered: θ Ct = max() where Z s x,1 is the thickness of the biocrust, and Z Ct is the thickness of the top crust.θ x was calculated from the surface humidity and the water retention of the crust layer using Eqs.( 25)-( 26) by Gong et al. (2016).

Calculation of litter input to soil and SOC transport in biocrust
The litterfall added to each soil layer (A 1 x,i ; Eq. 13) was linked to the mortality of roots, which was calculated following Asaeda and Karunaratne (2000).
where k mo is the optimal mortality rate at 20 • C; Q mo is the temperature sensitivity parameter (Asaeda and Karunaratne, 2000).Similarly, we attributed the C transport rate (A C m ) from M CA to M CH mainly to the mortality of autotrophic organisms.We assumed that most mortality of crust organisms occurred during abrupt changes in wetness, as microbial communities may adapt slow moisture changes or remain inactive during drought (e.g.Reed et al., 2012;Garcia-Pichel et al., 2013).Here, we introduced a water-content multiplier, f m (θ Ct ), to describe the impact of abrupt θ Ct changes on k m : where k mc is the optimal mortality rate at 20 • C; Q mo is the temperature sensitivity parameter (Asaeda and Karunaratne, 2000); θ Ct7 is the forward 7-day moving average of θ Ct .
C transport from top crust to subcrust was calculated as driven mainly by the sand deposition and burying of top crust SOM.Assuming the C content in top crust was homogeneous and the thickness Z Ct was near constant, the transport rate (k b ) was then proportional to the sand deposition rate: where ρ bulk is the bulk density of soil; k sand is the sand deposition rate in Yanchi area, which is a function of wind velocity (Li and Shirato, 2003).
J. Gong et al.: Modelling the diurnal and seasonal dynamics 2.5 Model validation and sensitivity analyses

Simulation setups
In the model simulations, soil depth was set to 67.5 cm to cover the rooting zone (Gong et al., 2016), including the crust layer (2.5 cm) and sandy subsoil (65 cm, stratified into 5 cm layers).Water contents measured at 70 cm depth was used as the lower boundary conditions for hydrological simulations (Jia et al., 2014).The calculation of soil temperature extended to 170 cm depth with the no-flow boundary, in regard to the probably strong heat exchange at the lower boundary of rooting zone (Gong et al., 2016).The zero-flow condition was set for the lower boundary of CO 2 and O 2 gases, whereas dissolved CO 2 was able to leach with seepage water.Based on presumed similarity of RLU structures, we assumed no-flux conditions for transport of water, heat, solvents and gases at the outer boundary.In the simulation, we assumed instant gas transport via top crust, whereas considered the CO 2 released by subcrust (R C s ) was subject to the dissolving-transport processes.In this work, we aggregated the C processes in subcrust with those in the soil profile.The initial ratio of M CA : M CH was set to 2 : 3. The C concentration of organic matter was set to 50 %.
The model was run with half-hourly meteorological variables including the incoming shortwave radiation, incoming longwave radiation, PAR, T a , relative humidity, wind speed and precipitation.Initial temperatures and soil moisture contents for each soil layer were initialized following the work of Gong et al. (2016).Surface CO 2 concentration was set to 400 ppm.The initial gaseous CO 2 concentration was set to increase linearly with depth (5 ppm cm −1 ).The initial CO 2 concentration in liquid form was then calculated based on Eqs. ( 4)-( 8).The initial content of mineral N content was set to 40 mg g −1 , which was within the range of the field observations.The two-dimensional transpiration of water, energy and gases along the soil profile was solved numerically using the predict-evaluate-correct-evaluate (PECE) method (Butcher, 2003).In order to avoid undesired numerical oscillations, the transport of water, energy and gases was calculated at 5 min substeps.

Model validation
First, we validated the modelling of soil temperature and moisture content for the F S site (Test 0).The simulated hourly soil temperature and moisture content at 10 cm depth were compared to the measured values for each collar.The validation was based on the same meteorological data as used by Gong et al. (2016), who validated the model in regard to the diurnal to seasonal dynamics of radiation balance, surface energy balance, soil temperature and moisture content at the EC site.
The validity of the modelled F S was then examined in three separate tests.In Test 1, modelled F S was validated for non-crusted soils.In this case, F T in Eq. ( 1) was the only term affecting F S (F B = 0 and F P = 0), and the crust influences on C-water exchanges were excluded.The biocrustrelated processes were considered in Test 2 and Test 3. Test 2 considered the dark respiration of biocrust (R Ct ) and set F B = R Ct and F P = 0. Test 3 considered all the flux components (F T , F P and F B ).In these tests, different values of root biomass were assigned to the model, regarding the different collar conditions (Table 1).In Tests 1-3, half-hourly F S was simulated and averaged to hourly, and compared to those measured from the collars C1-C3, respectively.Linear regressions were used to compare the modelled and measured values.The biases (ζ ) of the simulated values were calculated by subtracting the measured values from the modelled ones.Gap values in the measurements were omitted in the validation and the bias analyses.

Simulating componential CO 2 fluxes and their parameter sensitivities
Using the validated model, we simulated the temporal trends of C flux components (i.e.P Ct , R Ct , F P , F T , R a and R s ) in Test 4, in order to find out how the different flux components may have contributed to the total efflux (Table 3).The simulation used the same model setups and climatic variables as Test 3. It should be noted that although the model was built as an abstract for ecosystem-level processes, the simulation setups and validation were performed at a point level corresponding to respiration chambers.Therefore, understanding the uncertainty sourced from parameterization could be helpful for future development and applications.In Gong et al. (2016), we have studied the sensitivities of modelled soil temperature and moisture content to the variations in soil texture, water retention properties, vegetation parameters and plant-interspace heterogeneities.In this study, we also tested the sensitivity of F S and componential fluxes to the changes in a number of site-specific parameters (Table 4).These parameters included pH, nitrogen content, water-thermal conditions, root biomass, production rates and decomposition rates of litter, which are often key factors for regulating the soil C processes but likely to vary within and among ecosystems (see, e.g.Ma et al., 2011;Gong et al., 2016;Wang et al., 2015).Furthermore, we tested the model sensitivities to several newly defined parameters (i.e.n R , n p and f m ) in order to understand their effects on model uncertainties.F S and componential fluxes at interspace were simulated by varying the single parameter value by 10 or 20 %.The sensitivity of each tested flux was described by the difference (d F ) in the annual flux rate simulated using manipulated parameters, as compared to the rate simulated under no-change conditions.a R s + R a represents the total CO 2 production from soil respiration.R a is the total autotrophic respiration (R a = i R a i ; see Eq. 10) and R s is the total heterotrophic respiration (R s = i R s i ; see Eq. 12).b F C t − F P represents the net CO 2 exchanges of top crust; see Eqs. ( 17) and ( 24) for correspondent algorithms of the variables.For definitions of other fluxes, see Eq. ( 1) for F S , Eq. ( 3) for F T , Eq. ( 17) for F C t , Eq. ( 18) for P Ct and Eq. ( 24) for F P .
Table 4. Sensitivity of simulated F S and its componential fluxes to manipulations of parameter values.

Change of parameter
a Definitions of fluxes; see Table 3 and Sect

Comparing model sensitivities between plant cover and interspace
In order to study the effects of plant-interspace heterogeneity on soil CO 2 efflux, Test 5 simulated annual F S and componential fluxes at the plant cover and compared the values to interspace.The simulation setups were almost the same as those employed in Tests 1-3; the only exception was that the same initial values of SOC storages (650 gC m −2 ) and root biomass (200 g m −2 ) were used for under-canopy and interspace areas for comparison purposes.3 Results

Model validity
Compared to the EC site in previous work (Gong et al., 2016), the soil in this study was much coarser and the measured θ at 10 cm depth was constantly lower (Fig. 4), indicating the necessity of reparameterization and validation of the water-energy algorithms.Figure 4a shows the modelled hourly T s and θ at 10 cm depth with the mean values measured from the F S site during 2013.Based on the site-specific vegetation and soil texture parameters, our model explained 97 % of the variations in the measured hourly T s .The model underestimated the temperature mainly in summertime (i.e.days 150-250; Fig. 4a).The underestimation was more pronounced around the noontime in the diurnal cycle.As the water-content sensors may not accurately capture the moisture dynamics during the freezing period, only the simulations during the ice-free period were compared to measured data (Fig. 4b).During the ice-free period, the model explained 83 % of the variations in the measured mean water contents at 10 cm depth.The biases in the modelled temperature and moisture content were less than the spatial variations observed in this area (e.g.Wang et al., 2015).Therefore, our model could be able to reproduce the time series for the measured water-energy fluxes at the site.
Our model explained 87 and 83 % of the variations in the hourly F S measured on the non-crusted surface in the years 2013 and 2014 (Fig. 5a).The root mean square errors (RMSEs) were 0.43 and 0.29 µmol m −2 s −1 , respectively.The model mainly underestimated the daytime F S during the freezing periods.During the ice-free periods, the model mainly overestimated the efflux in early spring.The biases in modelling largely showed a diurnal pattern (Fig. 5b): F S was mainly underestimated in the afternoon hours (i.e. from 10:00 to 15:00 LT) but slightly overestimated in the afternoon and evening.At the daily level, our model explained 94 % of the variations in measured daily efflux during the 2-year period (Fig. 5c).
Compared to the non-crusted soil (C1), the simulated F S for crusted surfaces (C2 and C3) showed greater deviations from measured data.At the hourly scale, our model explained 75 % (year 2013) and 68 % (year 2014) of variations in measured F S from C2 (Fig. 6a), and 68 % (year 2013) and 61 % (year 2014) of variations in the F S measured from C3 (Fig. 6b).For the 2-year period, the RMSEs of the modelled hourly F S were 0.25 and 0.35 µmol m −2 s −1 for C2 and C3, respectively.The magnitudes of biases (|ζ |) were generally greater during the rainfall period (i.e. from the start of rainfall to 24 h after the end of rainfall) than the inter-rainfall period (Fig. 7).The simulated F S for C2 showed a similar diurnal pattern of biases as compared to C1, suggesting ineligible contributions from the biases in simulated subsoil emissions.Introducing photosynthesis and photodegradation of biocrust to the system (C3) led to greater overestimations in F S , and these were more obvious in the afternoon hours (i.e. from 12:00 to 18:00 LT) and during the wetting period.Nevertheless, at the daily scale, the model explained 91 % (C2; Fig. 5c) and 86 % (C3; Fig. 5d) of the variations in the measured F S during the 2-year period.There were no significant systematic deviations between the measured and the modelled daily values, as indicated by the regression slopes close to 1 and the intercepts close to 0 (Figs. 4 and 5).
The results above showed that the model was able to describe the seasonal variations of F S for both non-crusted and lichen-crusted soils.Moreover, the model captured the strong variability of hourly/daily F S in wetting-drying cycles.Compared to earlier statistical modelling by B. Wang et al. (2014) and W. Wang et al. (2014b), this model showed equal or improved accuracy.In this sense, we assume that our model has included the main mechanisms controlling the F S dynamics in the soil system and could be used for further analysis of componential C processes and their parameter sensitivities.

Modelled C flux components of F S
Test 4 showed that R s was the main contributor to the rootzone CO 2 production, which accounted for a major source of effluxes (F S ).Our measurements showed large diurnal  and seasonal variations in F S regardless the existence of crust covers (Figs. 5 and 6).Particularly, the F S dynamics depended strongly on rain events.Even at non-crusted soil (i.e.C1), F S dropped significantly from the pre-rainfall level even to near zero but rebounded rapidly and peaked after the rain stopped (Fig. 5a).This could relate to the mismatched trends of CO 2 production (R s + R a ) and emission (F T ) from the rooting zone with respect to the wetting-drying cycles (Fig. 8a).Compared to CO 2 production, the responses of F T to rainfall were generally lagged and smoothed (see examples in Fig. 8b-d), irrespective of the size of rain events.In the simulation, soil rewetting increased CO 2 production rapidly but depressed F T , which increased after rain ceased.In all the examples (Fig. 8b-d), F T exceeded R P within 48 h after the end of rain events.At the annual level, the total R P was larger during wetting period (i.e.rainfall days plus 1 day after rainfall) than during the rest of the days of the year (i.e. the drying period), whereas the total F T was greater during the drying period (Fig. 8e).
On the annual basis, CO 2 production (R s + R a ) and emission (F T ) from root-zone soil were mismatched (Table 3), and the former was more than 15 % greater than the latter.Such a gap was mainly due to the root uptake and transport of dissolved CO 2 (i.e.36 gC m −2 yr −1 ), whereas the loss of dissolved CO 2 via seepages or pore-mediated horizontal flows was limited (i.e.7.4 gC m −2 yr −1 ).The photosynthesis rate of top crust was 31.1 gC m −2 yr −1 at interspace.After rainfall, the C uptake by the top crust increased significantly and even turned the soil from net C source to sink during a few hours to 1 day (Figs.6, 8).However, at the annual scale, the C losses via respiration and photodegradation accounted for 90 % of the photosynthetic products, leading to a near-zero contribution of top crust to F S during the 2-year period (i.e.< 5 gC m −2 yr −1 ).
Analysis of parameter sensitivity showed that the modelled F S and the component fluxes were more sensitive to ±2 • C in T s or ±10 % in θ , as compared to the effects of ±10 or ±20 % in the other parameters (Table 4).Varying θ by 10 % produced greater impacts on the simulated R P and crust-related fluxes (i.e.P Ct , R Ct and F P ), as compared to changing T s by ±2 • C. Increasing θ by 10 % enhanced the simulated P Ct and F C t by 41 and 28 %, and doubled the net C sequestration (F C t -F P ) by the top crust.However, the contribution of such changes to annual F S was minor and accounted for only 2.0 % of the total efflux.Besides T s and θ , the simulated efflux was also sensitive to changes in root biomass (M R ).Manipulating root biomass by ±10 % changed the annual F T and F S by about 7 %, and such effects were 100 % greater than ±10 % changes in M tot in the soil.Adjustment of other parameters, e.g.n p (Eq. 31) and f m (Eq.38), had little influence on the modelled F S and the componential fluxes (Table 4).In addition, the model was robust to the adjustment of several crust-related parameters, i.e. k mc , M Ct and M CA : M CH .Hence, algorithms corresponding to those parameters could be simplified in future developments.

Modelled plant-interspace differences in C flux components
At either plant-covered or interspace area, R s was a major contributor to root-zone CO 2 produced, and F T dominated the total effluxes (Table 3).The C loss at interspace was 14 % faster than under canopy on an annual basis if root biomass and SOC were homogeneous at plant cover and interspace.The lower F S rate at plant cover is mainly attributed to the lower CO 2 production (R S +R a ) from subsoil.The C loss via seepage and root transport, which is the gap between subsoil CO 2 production and emission (F T ), was slightly higher under canopy (17 %) than at interspace (15 %).Compared to interspace, the photosynthesis of biocrust (P Ct ) was 34 % lower under canopy.This reduced the under-canopy F C t by 42 % in comparison with interspace.However, such a difference was largely offset by the reduced photodegradation rate under canopy, leading to limited plant-interspace differences in net sequestration by the top crust (i.e. by 1.4 gC m −2 yr −1 ).
We further compared the flux sensitivities at plant cover and interspace to the changes in three most effective parameters (i.e.T s , θ and M R ; see Table 5).For subsoil-mediated fluxes (i.e.F S , F T , R a , R s ), the sensitivity values differed by less than 2 % from plant cover to interspace.On the other hand, the sensitivities of crust-related fluxes (i.e.P Ct , F C t , F P ) showed greater differences between plant cover and interspace.Compared to interspace, F S and F T in plantcovered areas were more sensitive to T s changes but less sensitive to manipulations in θ .On the other hand, the plant cover reduced the sensitivity of CO 2 effluxes to changes in root biomass.P Ct , F C t and F P were generally more sensitive to warming and θ manipulations at plant cover than at interspace, except that plant cover decreased the sensitivity of F P to −10 % changes in θ .Nevertheless, their contribution to the sensitivity of F S was marginal due to the low flux rates of crusts.Our process-based model provided a tool to separate the multiple soil C processes and investigate their roles in regulating F S dynamics in dryland ecosystems.So far, efforts to quantify the soil C loss in terrestrial ecosystems have considered soil C efflux as a synonym of respired CO 2 .However, based on this work, precautions must be taken when extrapolating the F S responses from the chamber to ecosystem scale and from short-term to long-term periods.Processes other than autotrophic and heterotrophic respiration could significantly modify the F S responses to climatic variability.Our simulation highlighted decoupled CO 2 production and emission during the wetting-drying cycle, as regulated by the CO 2 transport in the soil profile.The simulated CO 2 production in the soil profile were much greater than effluxes during rain pulses (e.g.Fig. 7).This indicated that the low F S during rewetting was mainly due to the increase in CO 2 dissolution, instead of the reduced respiration rates by low O 2 supply (e.g.Fang and Moncrieff, 1999).This finding is further supported by the measurement of Maier et al. (2011), which showed that 40 % of the respired CO 2 could be stored temporally in the soil pore space after rainfall.The dissolved CO 2 was then released gradually with the evaporation of pore water, leading to lagged responses of efflux as compared to respiration.Given that a major fraction of CO 2 was produced during the wetting periods (Fig. 5e), such a lagging effect should be carefully examined when analysing the climatic sensitivity of F S .Our simulations showed that a considerable fraction of CO 2 produced could be removed by root uptake and leave the volume measured by the respiration chamber.Bloemen et al. (2016) showed that the CO 2 concentration in root xylems could be higher than in soil solutions.This implies that such a "missing source" might be even greater than the model estimation, although knowledge is still limited about the efficiency of the removal and the diffusion/release of CO 2 during the transport (Bloemen et al., 2016).The contributions of biocrusts as C sinks or sources have remained largely unknown (Castillo-Monroy et al., 2011).This is mainly due to the difficulty to separate the CO 2 exchanges of crust organisms from the background respiration (Castillo-Monroy et al., 2011;Sancho et al., 2016).As demonstrated in our work (Fig. 5b-d), the photosynthesis of top crust could be masked by background emission quickly (e.g.within 1 day) after rain events.The simulated F C t was 31 gC m −2 yr −1 at interspace.Considering a 30 % coverage of lichens over the sampling area (Feng et al., 2014), the interspace-level NPP was 9.3 gC m −2 yr −1 .This value was largely greater than the lab-based estimation for the site (Feng et al., 2014).However, it was in range of the values reported from several other dryland ecosystems (i.e.5.3-29 gC m −2 yr −1 ; Sancho et al., 2016).Our simulations also suggested that photodegradation might offset about 48 % of the CO 2 photosynthesized by biocrust.It could explain the much higher F S measured from the transparent chamber (C3) than from the opaque chamber (C2) during dry daytime periods (e.g.Fig. 9).It should be also noted that the litter from shrub canopy was not included in the measurement nor the modelling.Also, the interactions between photodegradation and biotic decay were not considered either.Hence, the contribution of photodegradation to soil C balance could be greater than our estimation at the ecosystem level (see, e.g.Gliksman et al., 2016).Although the contribution of surface exchanges was only marginal as compared to the annual CO 2 efflux, removing the biocrust processes would substantially reduce the model validity.For example, the goodness of fit (i.e.R 2 ) in Test 3 dropped from 0.65 to 0.45 for the 2-year period, if F C t and F P were neglected.Therefore, delineating the gas exchange of biocrust could be helpful in order to upscale the modelling of C balance from chamber to ecosystem level, where the distribution of the crust cover may vary from one site to another.

Plant-interspace differences in soil C fluxes
Clumped distributions of foliage and biomass are critical features for the adaptation and functioning of vegetation in arid and semiarid environments.Previous studies have mainly emphasized the shrub effects on ecohydrology (e.g.Gong et al., 2016) and enrichment of sediments and nutrients, known as the "resource island" effect (Reynolds et al., 1999;Rietkerk et al., 2004).Our simulations showed that the presence of shrub canopy also influenced soil C exchanges.The presence of shrub cover affected the C functioning of biocrust mainly through shading, which reduced photosynthesis more than respiration and photodegradation.Compared to interspace, the simulated annual F S was 13 % lower under canopy (Test 5).As we ruled out the differences in SOC and root biomass and limited the C-flux differences contributed by biocrusts between plant-covered and interspace areas, such a decrease in plant-cover F S was probably due to the cooling effect of canopy (Gong et al., 2016).This effect was close to the modelled responses of F S to ±2 • C in soil temperature or ±10 % in soil water content.As the root density and litter production rate are commonly larger under canopy than interspace (e.g.Zhang et al., 2008), the lower respiration rate under canopy tends to facilitate the accumulation of biomass and organic matter and feedback to the functioning of resource islands during prolonged periods.
Our simulation further indicated considerable differences in the C-flux sensitivities between areas with plant cover and those without.As the C processes and initial conditions were set to be homogeneous, those differences could mainly result from the different water-thermal conditions at plant cover and interspace.For example, the higher temperature sensitivities of F T , F S and F C t may relate to the cooling effect of canopy (see Gong et al., 2016), which may lead to a greater Q 10 value for respiration estimations (i.e.Eq. 27).Moreover, the slower decomposition in under-canopy soil could also lead to the lower sensitivities of F T and F S to changes in root biomass and SOM contents.On the other hand, water advection from interspace to plant cover, which may support over 30 % of water loss from under-canopy soils (Gong et al., 2016), could help to lower the F S sensitivity to water content changes at plant cover.The increased water-thermal sensitivities of C exchanges of biocrust could be explained by the less-stressful environment for curst organisms, e.g. higher moisture content but lower radiation and tempera-ture, although the photosynthesis of lichens (P Ct ) could be reduced by shading (Table 3).Such heterogeneity of C-flux sensitivities thus should be considered in future studies on the ecosystem-level responses to climate change and extreme climatic events.

Modelling uncertainties and future research needs
Our model showed its ability to describe the dynamics of soil temperature, moisture content and C effluxes measured for the studied semiarid ecosystem.Major uncertainty of the modelling, however, may refer to the concept of equifinality (Beven, 1993(Beven, , 2006)).The question of equifinality arises from the fact that the structures and mechanisms being modelled are based on insufficient information.Consequently, alternative models using different functions and parameter sets may fit equally well to observations; thus, the mechanisms quantified in modelling are difficult to be justified or falsified.Regarding this work, the modelling equifinality and uncertainty could relate to several aspects.
Firstly, the RLU was a statistical simplification to the target ecosystem at footprint scale (Gong et al., 2016) and may not fully capture the spatially explicit schemes of soil environment and biogeochemistry at ecosystem scale.For example, the model assumed Poisson probability of mutual shading (Bégué et al., 1994), and the probability of shading increased continuously with solar zenith (Gong et al., 2016).However, for explicit space-time, shading is binary.This possibly leads to the biases in the estimations of net radiation (Gong et al., 2016) and collar temperature around midday, which sequentially affected the simulated diurnal pattern of F S (see Fig. 3b).Moreover, field observations showed considerable spatial variations of soil temperature, water content and biogeochemistry (e.g.pH, litter quality and root biomass) within a distance of 3-5 m.Such variations could well exceed a magnitude of 10 % and even over 100 % (e.g.Zhang et al., 2008;Feng et al., 2013;Wang et al., 2015).Therefore, the variation of F S driven by the spatiality of soil factors could be greater than the responses to ±2 • C in soil temperature or ±10 % in soil water content.Therefore, future modelling may need to consider spatially explicit settings, in order to further minimize the gaps between model settings and reality.
Secondly, the high sensitivity of simulated F S to soil pH indicated that unconsidered processes of inorganic C could strongly affect the accuracy of modelling.Our modelling for CO 2 transport considered gaseous and liquid phases.However, the solid phase was not included, despite the high lime content (2300-5400 kg ha −1 ) in the soil (Feng et al., 2013;Wang et al., 2015).Based on soil samples of similar lime content (2700 kg ha −1 ), Buysee et al. (2013) showed that neglecting the inorganic C exchanges by solids may lead to underestimation of F S during the heating phase of the day but overestimation of F S during the cooling phase.This is very similar to the diurnal pattern of biases in simulated F S (Fig. 3d).Therefore, further improvement on the modelling may need to consider the solid phase as well.
Thirdly, the current model still lacked descriptions on growths of plant and soil microbes.Compared to many other ecosystems, drylands often feature a high root-shoot ratio (Jackson et al., 1996) but low SOC storage.Changes in plant physiology and growth can readily influence root metabolisms, and labile SOC pools hence modify F S dynamics (Wang et al., 2015).On the other hand, a large fluctuation of diurnal and seasonal temperature may drive the microbial communities to shift between those that are warm adapted and those that are cold adapted (Van Gestel et al., 2013), which could largely change soil respiration and its sensitivity to freeze-thaw cycles (Van Gestel et al., 2013;Liu et al., 2016).Both the biotic controls are mixed with the legacy effects of climatic variability over annual and interannual courses (Sala et al., 2012;Jia et al., 2016;Shen et al., 2016), and could affect the C-water simulations cumulatively through the feedbacks between biomass accumulation and soil biogeochemistry (Bradford et al., 2016).This may explain the decreasing trend of model validity from 2013 to 2014 (Figs. 3,4).Therefore, the dynamics of plant and microbial communities are required in future modelling in order to improve the F S simulations regarding interannual and long-term periods.
In addition, proper field data are still needed to support the future modelling work.The dataset used in our model validation mainly separated the influences of biocrusts from subsoil respiration.However, some processes like photodegradation and lateral CO 2 transport by root or water flows still require more support from observations.Also, respiration data from shrub-covered soil remain unavailable, as the settlement of soil collars and respiration chambers under canopy could easily interrupt the biophysiology of shrubs.The C functionality of crust organisms is especially sensitive to water content (Table 4).However, tracking the water content in the very thin layer of top crust can be very challenging using hourbased meteorological data.Nocturnal water inputs (e.g.dewfalls) are important for the metabolisms of crust organisms (e.g.Liu et al., 2006), but they are hard to be quantified precisely by EC measurement or models derived from EC data.Moreover, we presumed structural homogeneity for the crust layer and employed a constant regime for crust processes.In reality, there may not be clear boundaries between top crust and subcrust, and even top crust itself may contain significant variations in microstructure and communities even within 1 cm (Williams et al., 2012;Raanan et al., 2016).The C sequestration of biocrust can be strongly modified by microbial communities directly (Belnap et al., 2003;Pointing and Belnap, 2012;Feng et al., 2014;Maestre et al., 2015) and through other factors, e.g.surface albedo (Chamizo et al., 2012), dew falls (Liu et al., 2006) and soil pore formation (Williams et al., 2012;Felde et al., 2014).So far, many questions remain unanswered about the mechanisms that control the colonization, adaption and succession of microbial com-munities and the structure function of biocrust (Pointing and Belnap, 2012).Further knowledge on these mechanisms will be helpful to validate or falsify the modelled C functionality in response to climate change and extreme climatic events.

Conclusions
This work represents a first attempt to integrate the CO 2 production, transport and surface exchanges (e.g.biocrust photosynthesis, respiration and photodegradation) in F S modelling for dryland ecosystems with high plant-interspace heterogeneities.Our model reproduced the F S dynamics measured from non-crusted and lichen-crusted soil collars during 2013-2014, although introducing the gas exchanges of lichen crust decreased the model performance at the hourly scale.However, further model development may still be required on several aspects, e.g. by including (i) the spatially explicit schemes for surface conditions and soil biogeochemistry; (ii) influences of lime and solids on CO 2 transport; (iii) growth dynamics of plants; (iv) high-resolution dynamics of surface water-thermal conditions; and (v) the dynamics of microstructure and microbial communities of biocrusts.
Our model simulations highlighted that the transport processes of inorganic C and the metabolisms of biocrusts could strongly modify the CO 2 efflux, and these influences are closely linked to soil hydrology.Soil rewetting could enhance CO 2 dissolution and delay the emission of CO 2 produced from the root zone.In addition, an ineligible fraction of respired CO 2 could be removed via lateral flows and root uptake, and become "missing" from volumes under respiration chambers.The lichen-crusted soil could temporally shift from net CO 2 source to sink during rewetting, as driven by the photosynthesis of lichens and the restrained CO 2 emissions from subsoil, whereas after rain events the CO 2 exchanges of lichens could be easily masked by background emissions from subsoil.Based on our modelling, the annual NPP was 9.3 gC m −2 by the top crust at interspace.However, the net C sequestration by the top crust could be marginal if the photodegradation was accounted.Our modelling further showed different componential C fluxes and sensitivities between plant-covered and interspace areas.The presence of plant cover tended to decrease the root-zone CO 2 production and biocrust C sequestration but increase the temperature sensitivities of these fluxes.On the other hand, the sensitivities of root-zone emissions to water content were lower under canopy.This may be due to the advection water flows from the interspace to plant cover.To conclude, the complexity and plant-interspace heterogeneities of soil C processes should be carefully considered when extrapolating findings from chamber to ecosystem scales, in order to predict the ecosystem responses to climate change and extreme climatic events.Our model can serve as a useful tool to simulate the soil CO 2 efflux dynamics in dryland ecosystems.

Figure 2 .
Figure 2. Conceptual framework of process-based modelling.Solid arrows represent flows of masses and dashed arrows represent flows of information.

Figure 3 .
Figure 3. Measured and fitted bulk respiration (a) and photosynthesis (b) of the lichen top crust as functions of temperature and water content.
. 2.5.3.b Value represents the percentage (%) of change (d F ) in correspondent C flux with manipulated parameter value, as compared to the no-change condition.A positive value represents the percentage of increase in the simulated flux, whereas a negative value represents the percentage of decrease.c The change in simulated C flux was smaller than 1 %.

Figure 4 .
Figure 4. Measured and modelled soil temperature (a) and soil moisture content (b) at 10 cm depth for the F S site and as compared to the EC site in the year 2013 by Gong et al. (2016).

Figure 5 .
Figure 5. Measured and modelled hourly F S for non-crusted soil (a), the temporal pattern of the bias of simulated hourly F S (b) and the comparison of measured and modelled daily F S (c) during 2013-2014.

Figure 6 .
Figure 6.Measured and modelled F S of lichen-crusted soils for opaque (a, c) and transparent chambers (b, d) at hourly (a, b) and daily (c, d) scales during 2013-2014.

Figure 7 .
Figure 7.Diurnal patterns of biases (ζ ) in the simulated hourly F S for lichen-crusted soils using opaque (a) and transparent chambers (b), and the cumulative probability of the biases during wetting and drying periods (c) during 2013-2014.The wetting period included the rainy days and a 1-day period after each rainfall.The drying period included the rest of the years other than the wetting period.
roles of componential C processes in regulating soil CO 2 efflux

Figure 8 .
Figure 8. Simulated component CO 2 exchanges by biocrust and root-zone soil (a), the simulated CO 2 fluxes before and after example rain events of 2.3 mm (b), 7.6 mm (c) and 12.8 mm (d) sizes, and the comparison of F T and R P during wetting and drying periods during 2013-2014 (e).The wetting period included the rainy days and a 1-day period after each rainfall.The drying period included the rest of the years other than the wetting period.

Figure 9 .
Figure 9.Comparison of the measured F S from lichen-crusted surfaces using opaque and transparent chambers during a dry period (days 83-103) in spring 2013.

Table 1 .
Configuration of soil collars used in this study.
*The values were calculated from the measured hourly F S data excluding data gaps.

Table 5 .
Plant-interspace differences in the sensitivities of C fluxes to changes in soil temperature (T s ), water content (θ) and root biomass (M R ).
Table3 and Sect.2.5.3.bValuesshows the plant-interspace difference in parameter sensitivities by value (outside bracket, d Fp − d Fi ) and by percentage (inside bracket, 100 × (|d Fp | − |d Fi |)/|d Fi |), where d Fp and d Fi are parameter sensitivities (d F ; for definitions, see Table4and Sect.2.5.3 ) for plant-covered and interspace areas, respectively.A positive percentage (inside bracket) indicates a greater sensitivity (|d F |) of the flux at the plant cover than at interspace, whereas a negative value indicates a lower sensitivity.For definitions of fluxes and sensitivities, see Table3 and Sect.2.5.3.cThedifference in sensitivity is smaller than 0.1 % by value.