How vadose zone mass and energy transfer physics affects the ecohydrological dynamics of a Tibetan meadow?

The vadose zone is a sensitive region to environmental changes and exerts a crucial control in ecosystem functioning. While the way in representing the underlying process of vadose zone differs among 15 models, the effect of such differences on ecosystem functioning is seldomly reported. Here, the detailed vadose zone process model STEMMUS was coupled with the ecohydrological model T&C to investigate the role of solving influential physical processes, considering different soil water and heat transfer parameterizations including frozen soils. We tested model performance with the aid of a comprehensive observation dataset collected at a typical meadow ecosystem on the Tibetan Plateau. Results indicated that: 20 i) explicitly considering the frozen soil process significantly improved the soil moisture/temperature (SM/ST) profile simulations and facilitated our understanding of the water transfer processes within the soil-plantatmosphere continuum; ii) the difference among various complexity of vadose zone physics have an impact on the vegetation dynamics mainly at the beginning of the growing season; iii) models with different vadose zone physics can predict similar interannual vegetation dynamics, and energy, water and carbon exchanges 25 at the land-surface. This research highlights the role of vadose zone models and their underlying physics, in ecosystem functioning and can guide the development and applications of future earth system models. 30 https://doi.org/10.5194/tc-2020-88 Preprint. Discussion started: 12 May 2020 c © Author(s) 2020. CC BY 4.0 License.


Introduction
Understanding how ecosystem functioning interacts with changing environmental conditions is a crucial yet challenging problem of earth system research. Various types of models, including land surface models, terrestrial biosphere models, ecohydrology models, and hydrological models, have been widely utilized to 35 enhance our knowledge in terms of land surface and hydrological processes including the role of vegetation (Fatichi et al., 2016a;Fisher et al., 2014). A number of uncertainties are originated from different model structures. A more detailed knowledge of the effect of using a given model formulation can help toward making better projections of land surface dynamics, also in response to the call for joint efforts for systematic model developments (Clark et al., 2015;Yu et al., 2016). With emphasis on enhancing the underlying physics, 40 most of the models have already adopted various solutions and parameterizations of land surface processes (e.g., Noah MP, CLM5, T&C) (Fatichi et al., 2012a, b;Niu et al., 2011;Lawrence et al., 2019), which facilitate the appropriate descriptions of different physical processes in various ecosystems. However, in these models, the water and heat transfer process in the vadose zone remains independent and uncoupled, as they often adopt simplified approaches to water and heat transfer in the subsurface. Such physical parameterizations of 45 vadose zone might result in unsatisfactory simulations or physical interpretations, especially when water and heat are tightly coupled as for instance in freezing soils (Hansson et al., 2004). In this regard, researchers have stressed the necessity to simultaneously consider the water and heat transfer process in dry/cold seasons (Bittelli et al., 2008;Scanlon and Milly, 1994;Yu et al., 2016;Yu et al., 2018;Zeng et al., 2009a;Zeng et al., 2009b).

50
With the largest area of high-altitude permafrost and seasonally frozen ground, Tibetan Plateau is recognized as one of the most sensitive regions for climate change (Cheng and Wu, 2007;Liu and Chen, 2000).
Monitoring and projecting the dynamics of hydrothermal and ecohydrological states and their responses to climate change in the Tibetan Plateau is considerably important to help shedding light on future ecosystem responses. Considerable land-surface and vegetation changes have been reported, e.g., degradation of 55 permafrost and changes of frozen ground (Cheng and Wu, 2007), advancing vegetation leaf onset dates (Zhang et al., 2013), and enhanced vegetation activity at start of growing season (Qin et al., 2016). However, there are divergences with regard to the expected ecosystem modifications across the Tibetan Plateau (Qin et al., 2016;Wang et al., 2018;Zhao et al., 2010). It is thus fundamental to have in situ multicomponent measurement networks (including meteorology, soil moisture/temperature, surface energy fluxes, carbon 60 fluxes) to understand the environmental controls (Hao et al., 2011;Wang et al., 2018;Wang et al., 2017;Zhao et al., 2010), validate terrestrial biosphere models and remote sensing products (He et al., 2014;Mwangi et al., 2020;Niu et al., 2016;Su et al., 2013;Tian et al., 2017), and extrapolate results via model-data-fusion methods to larger scales to better characterize land surface processes and ecosystem dynamics of the Tibetan Plateau (He et al., 2014;Zeng et al., 2016;Zhuang et al., 2020).

65
In this study, we tested the consequences of considering coupled water and heat transfer processes on landsurface fluxes and ecosystem dynamics in the extreme environmental conditions of the Tibetan plateau https://doi.org/10.5194/tc-2020-88 Preprint. Discussion started: 12 May 2020 c Author(s) 2020. CC BY 4.0 License.
relying on state-of-the-art land-surface and ecohydrological modeling confronted with multiple field observations. The limited knowledge of including or not complex vadose zone processes in such environment frames the scope here. Specifically, the driving questions of the research are: i) How different complexity in 70 representing frozen soil and coupled water and heat physics is affecting the simulated ecohydrological dynamics of a Tibetan plateau meadow? ii) How does model complexity affect our interpretation of mass, energy, and carbon fluxes at the ecosystem scale? Answering these questions is important to evaluate the adequacy of models to answer questions related ecosystem changes across the Tibetan Plateau.
In order to achieve the aforementioned goals, the detailed soil mass and energy transfer process developed 75 in the STEMMUS model (Zeng et al., 2011a, b;Zeng and Su, 2013) was incorporated into the ecohydrology model Tethys-Chloris (T&C) (Fatichi et al., 2012a, b). The frozen soil physics was explicitly taken into account and soil water and heat transfer are fully coupled to further facilitating the model's capability in dealing with complex vadose zone processes.

Experimental site
The Maqu soil moisture and soil temperature (SMST) monitoring network (Dente et al., 2012;Su et al., 2011;Su et al., 2013;Zeng et al., 2016) is situated on the north-eastern fringe of the Tibetan Plateau. The monitoring network covers an area of approximately 40 km×80 km (33°30'-34°15'N, 101°38'-102°45'E) with the elevation varying from 3200 m to 4200 m above the sea level (a.s.l.). The climate can be 85 characterized by wet rainy summers and cold dry winters. The mean annual air temperature (MAT) is 1.2℃ with about -10.0℃ and 11.7℃ for the coldest month (January) and warmest month (July), respectively. The alpine meadows (e.g., Cyperaceae and Gramineae) dominate in this region with the height of about 5 cm during the wintertime and 15 cm during the summertime. The general soil types are categorized as sandy loam, silt loam with a maximum of 18.3 % organic matter for the upper soil layers (Dente et al., 2012;Zhao 90 et al., 2018;Zheng et al., 2015a;Zheng et al., 2015b). The groundwater level of the grassland area fluctuates from about 8.5 m to 12.0 m below the ground.
At Maqu site, SMST profiles (5 cm, 10 cm, 20 cm, 40 cm, and 80 cm) are automatically measured by 5 TM ECH2O probes (METER Group, Inc., USA) at a 15-min interval. The meteorological forcing (including wind speed/direction, air temperature and relative humidity at five heights above ground) is recorded by a 20 m 95 Planetary Boundary Layer (PBL) tower system. An eddy-covariance system (EC150, Campbell Scientific, Inc., USA) was installed for monitoring the dynamics of the turbulent heat fluxes and carbon fluxes.
MCD15A3H provides estimation of 8-day composites of LAI and FPAR, while MOD17A2H an 8-day composite of Gross Primary Production (GPP). Both MODIS products are at a resolution of 500m.

Land surface carbon fluxes
Starting from the raw NEE (Net Ecosystem Exchange) and ancillary meteorological data (friction velocity * , global radiation , soil temperature , air temperature , and vapor pressure deficit ), we employed the REddyProc package (Reichstein et al., 2005;Wutzler et al., 2018) as post-processing tool to obtain the time series of NEE, GPP (Gross Primary Production) and ecosystem respiration dynamics.

110
Three different techniques, * filtering, gap filling, and flux partitioning, were adopted in REddyProc package. The periods with low turbulent mixing is firstly determined and filtered for quality control ( * filtering, (Papale et al., 2006)). Then, the marginal distribution sampling (MDS) algorithm was used as the gap filling method to replace the missing data (Reichstein et al., 2005). Finally, NEE was separated into GPP and by night-time based and day-time based approaches (Lasslop et al., 2010).

Precipitation, evapotranspiration, and frost front
The observed surface water conditions over the entire study period, including the precipitation and cumulative evapotranspiration (which is obtained by summing up the hourly latent heat flux measured by EC system), are shown in Fig. 1a. Both ET and precipitation are low until the end of the freezing period (see Fig.   1b), during this early period the daily average ET is 0.15 mm/d. During the growing season, the cumulative 120 precipitation increases and ET follows with a lower magnitude. The average daily ET for the entire observation period is 1.45 mm/d.

130
The Tethys-Chloris model (T&C) (Fatichi et al., 2012a) simulates the coupled dynamics of energy, water, and vegetation and has been successfully applied to a very large spectrum of ecosystems and environmental conditions as summarized elsewhere (Fatichi and Ivanov, 2014;Fatichi and Pappas, 2017;Fatichi et al., 2016b;Mastrotheodoros et al., 2017;Pappas et al., 2016). The model simulates the energy, water, and carbon exchanges between the land surface and the atmospheric surface layer accounting for aerodynamic, references and some key elements for this article are discussed in the following.

Overview of STEMMUS
STEMMUS (Simultaneous Transfer of Energy, Mass and Momentum in Unsaturated Soil) model solves soil water movement, soil air flow, and soil heat flow balance equations simultaneously in one timestep (Zeng et al. 2011a,b;Zeng and Su, 2013). The Richards' equation with modifications made by Milly (1982) was 145 utilized to mimic the coupled soil mass and energy transfer process. The dry air is considered to be an independent phase in the soil. The vapor diffusion, advection and dispersion are all taken into account as the water vapor transport mechanism. In addition to the soil moisture and temperature gradient, the atmospheric pressure gradient acts as the third driving force for soil water, vapor and heat flow. Root water uptake process is regarded as the sink term of soil water and heat balance equations, building up the linkage between soil 150 and atmosphere (Yu et al., 2016).

Difference between STEMMUS and T&C models in mass and energy transfer process
While T&C model specializes in dealing with the interaction between vegetation and the hydrological system, it simplified the soil water and heat transfer process in the hydrology component, e.g., ignored water vapor flow, dry air flow, and, in the original version, does not contain freezing/thawing process, as water is always 155 in liquid phase regardless of sub-zero soil temperatures. To extend the application of T&C model over frozen soil, a freeze-thaw module has been incorporated for this study as described below. Furthermore, while STEMMUS model can well reproduce the soil water and heat transfer process, it lacks detailed description of land-surface processes and of the vegetation-hydrology feedback mechanisms. To take advantage of the strengths of both models, we coupled STEMMUS with the land-surface and vegetation components of T&C 160 model to better describe the soil-plant-atmosphere continuum.

1) Mass transfer process
The 1-D Richards equation, which describes the water flow under gravity and capillary forces in isothermal conditions, is solved in T&C for variably saturated soils.
where (m 3 m -3 ) is the volumetric water content; q (kg m -2 s -1 ) is the water flux; z (m) is the vertical direction 165 coordinate; S (s -1 ) is the sink term for transpiration, evaporation and lateral transfer fluxes. ρL (kg m −3 ) is the liquid water density; K (m s -1 ) is the soil hydraulic conductivity; (m) is the soil water potential; t (s) is the time.
In T&C, the nonlinear partial differential equation is solved using a finite volume approach with the method of lines (MOL) (Lee et al., 2004). MOL discretizes the spatial domain and reduces the partial differential 170 equation to a series of ordinary differential equations in terms of time, which can be expressed as where ( where ρV and ρi (kg m −3 ) are the density of water vapor and ice, respectively; θL , θV and θice (m 3 m −3 ) are the soil liquid, vapor and ice volumetric water content, respectively; ℎ , , and (kg m -2 s -1 ) are the soil 180 liquid water flow driven by the gradient of soil matric potential , temperature , and air pressure , respectively. ℎ , , and (kg m -2 s -1 ) are the soil water vapor fluxes driven by the gradient of soil matric potential , temperature , and air pressure , respectively. T (°C) is the soil temperature; and Pg (Pa) is the mixed pore-air pressure.
(kg m -2 s -2 ) is the specific weight of water. DTD (kg m -1 s -1 °C -1 ) is the transport coefficient of the adsorbed liquid flow due to temperature gradient; DVh (kg m -2 s -1 ) is the isothermal 185 vapor conductivity; and DVT (kg m -1 s -1 °C -1 ) is the thermal vapor diffusion coefficient; DVa is the advective vapor transfer coefficient.

2) Energy transfer process
The heat conservation equation used in the original T&C neglects the coupling of water and heat transfer physics and only the heat conduction component is considered, which can be expressed as below where Csoil (J kg −1 °C −1 ) is the specific heat capacities of bulk soil; λeff (W m −1 °C −1 ) is the effective thermal conductivity of the soil. When soil undergoes freezing/thawing process, the latent heat due to water phase change becomes important but is not included in the original T&C model. heat of freezing/thawing (− ) and a source term associated with the exothermic process of wetting of a porous medium (integral heat of wetting) (− ).
where ρs (kg m −3 ) is the soil solids density; θs is the volumetric fraction of solids in the soil; Cs, CL, CV, Ca and Ci (J kg −1 °C −1 ) are the specific heat capacities of soil solids, liquid, water vapor, dry air and ice, respectively; Tr (°C) is the arbitrary reference temperature; L0 (J kg −1 ) is the latent heat of vaporization of 200 water at the reference temperature; Lf (J kg −1 ) is the latent heat of fusion; W (J kg −1 ) is the differential heat of wetting (expressed by Edlefsen and Anderson (1943) as the amount of heat released when a small amount of free water is added to the soil matrix). qL, qV, and qa (kg m -2 s -1 ) are the liquid, vapor water flux and air flow, respectively. Additional details and the air flow balance equation for solving the coupled water and heat equations can be found in Zeng et al. (2011a, b) and Zeng and Su (2013). 205

T&C model with freezing/thawing process
In the T&C version modified to explicitly account for freezing/thawing processes, the heat conservation equation is written as: where the latent heat associated with the freezing/thawing process is explicitly considered and ice water content θi is a prognostic variable, which is simulated along with liquid water content for each soil layer.

210
Specifically, when Eq. (7) is rewritten in terms of an apparent specific heat capacity Capp (Gouttevin et al., 2012;Hansson et al., 2004) it can be solved equivalently to Eq. (4): where Capp can be computed knowing the temperature T, latent heat of fusion Lf and the differential (specific) water capacity dθ/dψ at a given water content (Hansson et al., 2004): The effective thermal conductivity λeff (W m −1 °C −1 ) and the specific soil heat capacity Csoil (J kg −1 °C −1 ) are 215 computed accounting for solid particles, water, and ice content (Farouki, 1981;Johansen, 1975;Yu et al., 2018;Lawrence et al., 2018). The soil freezing characteristic curve providing the maximum liquid water content at a given temperature is computed following Dall'Amico et al. (2011) and it can be combined with various soil hydraulic parameterization including van-Genuchten, Clapp and Hornberger and Saxton and Rawls (Fuchs et al., 1978;Yu et al., 2018). 220 https://doi.org/10.5194/tc-2020-88 Preprint. Discussion started: 12 May 2020 c Author(s) 2020. CC BY 4.0 License.
Finally, saturated hydraulic conductivity is corrected in presence of ice content (e.g., (Hansson et al., 2004;Yu et al., 2018)). Note, that beyond latent heat associated with phase change and changes in thermal and hydraulic parameters because of ice presence, all the other soil physics processes described by STEMMUS are neglected and heat and water fluxes are still uncoupled in this version of T&C.

225
The current coupling procedure between STEMMUS and the original T&C is based on a sequential coupling via exchanging mutual information within one timestep (see Figure 2).

250
All three numerical experiments shared the same soil and vegetation parameter settings to accommodate the conditions of a Tibetan meadow. The total depth of soil column was set as 3m and divided into 18 layers with a finer discretization in the upper soil layers (1-5cm) than that in the lower soil layers (10-50cm). Soil samples were collected and transported to the laboratory to determine the soil hydrothermal properties (see Zhao et al. (2018) for detail). The average soil texture and fitted Van Genuchten parameters at three soil layers were listed in supplement Table S1. Vegetation parameters were obtained on the basis of physical constraints, literature, and expert knowledge (see a summary of the adopted vegetation parameters in the supplement Table S2).

260
The 5-day moving average dynamics of the net incoming radiation (Rn), latent heat (LE) and sensible heat (H) fluxes measured and simulated by the unCPLD model, unCPLD-FT and CPLD model for the study period are presented ( Figure 3). The seasonality and magnitude of surface fluxes can be captured across seasons. A good match between observed and simulated Rn and LE was identified during the whole period, with isolated observable discrepancies (Fig. 3a &c). These mismatches of Rn can be partly attributed to the 265 uncertainties of observed winter precipitation events and the following snow cover dynamics, which might not be well captured in the models, because the true winter precipitation is difficult to observe. For the sensible heat flux simulations, all three models can reproduce the seasonal dynamics. However, an overestimation of the 5-day average values was observed in several periods. Given the good correspondence between observations and simulations of net radiation and latent heat, this discrepancy might be a model 270 shortcoming but it can be also generated by the lack of energy balance closure in the flux tower data (see Sect. 4.1). Compared with unCPLD and unCPLD-FT simulations, the overestimation was however reduced by the CPLD model simulations.
The correlation between observed and simulated surface heat fluxes with unCPLD, unCPLD-FT, and CPLD model is shown in Fig. 4. Noticeably all the unCPLD/CPLD model scenarios, with different water and heat 275 transfer physics, exhibit nearly identical statistical performance of surface fluxes simulations (Fig. 4). There is an overestimation of H reproduced by three model simulations. The CPLD model presented less overestimation of H compared to unCPLD models. The overall performance of the model in terms of turbulent flux simulations could be regarded as acceptable given current uncertainties in winter precipitation and flux-tower observations in such a challenging environment, even though discrepancies exist during 280 certain periods (Fig. 3).

Soil moisture simulations
The capability of the three models to reproduce the temporal dynamics of soil moisture is illustrated in Figure   5. By explicitly considering soil ice content, the unCPLD-FT and CPLD model capture well the response of soil moisture dynamics to the freezing/thawing process. While the unCPLD model lacks such capability and 285 maintains a higher soil moisture throughout the winter period, which then reflect in slightly lower water contents in the growing season. For all three models, the consistency between the measured and model simulated soil water content at five soil layers is satisfactory during the growing season, indicating the models' capability in portraying the effect of precipitation and root water uptake on the soil moisture conditions. https://doi.org/10.5194/tc-2020-88 Preprint. Discussion started: 12 May 2020 c Author(s) 2020. CC BY 4.0 License.

290
Five layers of soil temperature measurements were employed to test the performance of model in reproducing the soil thermal regime profiles (Fig. 6). During the growing period, all three models can well capture the dynamics of soil temperature at various depth with fluctuating atmospheric forcing. In this period, there is no significant difference among the three models with regards to the magnitude and temporal dynamics of soil temperature. During the freezing period, a general underestimation of soil temperature and 295 overestimation of its diurnal fluctuations were found at shallower soil layers, which may indicate that there is some thermal buffering effect in reality not fully considered in the models. Compared with unCPLD-FT and CPLD models, the unCPLD model simulations have stronger diurnal fluctuations of soil temperature with an underestimation of temperature at the beginning of the freezing period and a considerable overestimation during the thawing phase. This results in an earlier date passing the 0°C threshold than in the 300 unCPLD-FT simulations.

Soil ice content and water flux
The simulated time-series of soil ice content and water flux from both unCPLD and CPLD model simulations for soil layers below 2 cm are presented in Figure 7. As soil ice content measurements were not available, the freezing front propagation inferred from the soil temperature measurements was employed to qualitatively 305 testify the model performance. A deeper presence of soil ice content was reproduced by the unCPLD-FT model, as indicated by a deeper freezing front propagation than the in-situ measurements. CPLD model presented a relative good match of soil freezing dynamics as it is physically constrained by the interdependence of liquid, ice, vapor, air components in soil pores. The phenomenon that a certain amount of liquid water flux moves upwards along with the freezing front can be clearly noticed from the unCPLD-FT 310 and CPLD model simulations. As the soil matric potential changes sharply during the water phase change period, a certain amount of water fluxes will be forced towards the phase changing region, a phenomenon known as cryosuction. Such a phenomenon has already been demonstrated from theoretical and experimental perspectives by many researchers (Hansson et al., 2004;Watanabe et al., 2011;Yu et al., 2018). This is of course absent from the unCPLD model simulations (Fig. 7c). Nevertheless, the precipitation induced 315 downward water flux can be observed in both models.

Simulations of land surface carbon fluxes
The eddy covariance derived and remote sensing (MODIS) observations of vegetation dynamics are compared with the model simulation in Fig. 8. When compared with in situ flux-tower observations, a slightly earlier growth and considerably earlier senescence of grassland with lower photosynthesis was inferred from 320 MODIS GPP product (Fig. 8a).  Table S2). While the EC system used in this experimental site observes an annual GPP and NEE as 1132.52 and -293.24 g C m -2 yr -1 . Both the GPP and NEE measured fluxes are significantly larger than previous estimates of carbon exchange for such a type of ecosystem (and representative of much more productive ecosystems) and are unlikely to be correct in magnitude. The ecosystem respiration ( ), indicating the respiration of activity of all living organisms in an ecosystem is 340 shown in Fig. 8d. The performance of all three model simulations in reproducing dynamics can be characterized as an overall good match with regards to the magnitude and seasonal dynamics of , which further suggest the discrepancy in observed/simulated GPP is the driver of the discrepancy of NEE.
The difference in the soil liquid water/temperature profile simulations between the CPLD and unCPLD models (as shown in Fig. 5 & 6) resulted in difference in simulated vegetation dynamics, especially 345 concerning the leaf onset date, which is affected by integrated winter soil temperatures. The unCPLD simulations have a slightly lower vegetation activity compared to the CPLD model simulations either because of a delayed leaf onset (unCPLD-FT) or because of a slightly enhanced water-stress (unCPLD) induced by the different soil-moisture dynamics during the winter season. Indeed, there is a slight lower root zone moisture produced by the unCPLD model (Fig. 5), which affects the plant photosynthesis and growth, thus 350 the vegetation dynamics in T&C. The unCPLD-FT model has a delay in the vegetation onset date, due to its prolonged freezing conditions as derived from soil temperature simulations than the other simulations.

Surface energy balance closure
The energy balance closure problem, usually identified because the sum of latent (LE) and sensible (H) heat 355 fluxes is less than the available energy (Rn-G0), is quite common in eddy covariance measurements (Leuning et al., 2012;Wilson et al., 2002). The energy imbalance of EC measurements is particularly significant at the sites over the Tibetan Plateau (Tanaka et al., 2003;Yang et al., 2004;Zheng et al., 2014). Figure 9 presents https://doi.org/10.5194/tc-2020-88 Preprint. Discussion started: 12 May 2020 c Author(s) 2020. CC BY 4.0 License.
the energy balance imbalance of hourly LE and H by the eddy covariance measurements, observed Rn by the four-component radiation measurements, and the estimated ground heat flux (G0) by CPLD model. The sum 360 of measured LE and H was significantly less than Rn, with the slope of 0.59 (Fig. 9a). Usually, the measurements of radiation are reliable (Yang et al., 2004). If we assume that the turbulence fluxes (LE, H) measurements were accurate, then the rest of energy (around 41% of Rn) should be theoretically consumed by ground heat flux G0, which is clearly impossible. When compared to the available energy (Rn-G0), the slope was increased to 0.70 (Fig. 9b). Table 1 further demonstrated that the energy imbalance problem was 365 significant across all seasons. The seasonal variation of energy closure ratio (ECR) can be identified for the case LE+H versus Rn-G0, similar to the research of Tanaka et al. (2003), i.e., a good energy closure during the pre-monsoon periods while a degraded one during the summer monsoon periods.
These problems are clearly suggesting that care should be taken in the model to data comparison, but they are not affecting the comparison among models with different complexity as we did not force any parameter 370 calibration or data-fitting procedure, but simply rely on physical constraints, literature, and expert knowledge to assign model parameters.

Effects on water budget components
The effect of different model scenarios on soil water budget components is illustrated in Fig. 10. surface evaporation process (evaporation from soil surface and intercepted water). ii) water storage amount 395 in the vadose zone increased while the bottom leakage decreased. We attribute this to the way ice content is simulated in the CPLD simulation, and also to the temperature dependence of soil hydraulic conductivity.
Taking into account the fully coupled water and heat physics modify the temporal dynamics of ice formation and thawing in the soil and activates temperature effects on water flow (i.e., low soil temperature will slow down water movement). That implies that soil water flow toward and at the bottom soil layer is retarded 400 when coupled water and heat physics is considered (as reflected by less leakage water flux for CPLD simulations).

The potential influential pathways of different mass/heat transfer processes
Given the same atmospheric forcing and the same model structure to represent land-surface exchanges and vegetation dynamics, how water/heat transfer processes are represented in the soil generates differences in 405 SM and ST vertical profiles. From the perspective of energy and carbon fluxes, the convective heat flux and explicit frozen soil physics are taken into account in the CPLD model while they are not considered in the unCPLD models. The liquid water flux induced convective heat flux is mostly relevant during the frozen period. As it has been observed, a certain amount of liquid water/vapor flux moving toward the freezing front ( Fig. 7), resulting in a convective heat toward the front. Such amount of heat and mostly the heat release by 410 freezing and consumed by the melting processes slows down the freezing/thawing process and decreases the diurnal and seasonal temperature fluctuations (Fig. 5). Different soil thermal profiles have consequences on the vegetation dynamic process (Fig. 9), mainly by modifying the temperature profile in the soil, which affects the beginning of the growing season and the subsequent simulated photosynthesis and growth processes. From the perspective of water fluxes, it is during the frozen period that water and heat transfer 415 process are tightly coupled. Both the explicit consideration of soil ice and coupled water and heat physics can affect the vadose zone water flow via altering the hydraulic conductivity. This is testified by the fact that even the unCPLD-FT simulation accounting for soil-freezing in a simplified way in comparison to STEMMUS (e.g., the CPLD simulation) cannot recover the exact same dynamics of ice content (Fig. 7), which impacts leaf onset and to a less extent hydrological fluxes. However, in the rest of the year the 420 simplified solution of vadose zone physics of T&C leads to very similar results as the coupled one, suggesting that most of the additional physics does not modify substantially the ecohydrological response.

Conclusion
The

455
Competing interests. The authors declare that they have no conflict of interest.