The Role of Vadose Zone Physics in the Ecohydrological Response of a Tibetan Meadow to Freeze-Thaw Cycles

The vadose zone is a sensitive zone to environmental changes and exerts a crucial control in ecosystem functioning, and even more so in cold regions considering the rapid change of seasonally frozen ground under climate warming. While the way in representing the underlying physical process of vadose zone differs among models, the effect of such differences on ecosystem functioning and its ecohydrological 15 response to freeze-thaw cycles is seldom reported. Here, the detailed vadose zone process model STEMMUS was coupled with the ecohydrological model T&C to investigate the role of influential physical processes during freeze-thaw cycles. The physical representation is increased from using T&C without, and with explicit consideration of the impact of soil ice content on energy and water transfer properties, to T&C coupling with STEMMUS enabling the simultaneous mass and energy transfer in the soil system (liquid, 20 vapor, ice). 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: i) explicitly considering the frozen soil process significantly improved the soil moisture/temperature profile simulations and facilitated our understanding of the water transfer processes within the soil-plant-atmosphere continuum; ii) the difference among various representations of vadose zone physics have an impact on the vegetation dynamics mainly at 25 the beginning of the growing season; iii) models with different vadose zone physics can predict similar interannual vegetation dynamics, energy, water, and carbon exchanges at the land surface. This research highlights the important role of vadose zone physics for ecosystem functioning in cold regions and can support the development and application of future Earth system models.


Introduction
Recent climatic changes have accelerated the dynamics of frozen soils in cold regions, as for instance favoring permafrost thawing and degradation (Cheng and Wu, 2007;Hinzman et al., 2013;Peng et al., 2017; 35 Yao et al., 2019;Zhao et al., 2019). In consequence of these changes, vegetation cover and phenology, land surface water and energy balances, subsurface soil hydrothermal regimes, water flow pathways were reported to be affected (Wang et al., 2012;Schuur et al., 2015;Walvoord and Kurylyk, 2016;Gao et al., 2018;Campbell and Laudon, 2019;Yu et al., 2020). Understanding how an ecosystem interacts with changing environmental conditions is a crucial yet challenging problem of Earth system research for high 40 latitude/altitude regions and deserves further attention.
Land surface models, terrestrial biosphere models, ecohydrology models, and hydrological models have been widely utilized to enhance our knowledge in terms of land surface processes, ecohydrological processes (Fatichi and Ivanov, 2014;Fatichi et al., 2016a), and freezing/thawing process (Ekici et al., 2014;Cuntz and Haverd, 2018;Wang and Yang, 2018;Druel et al., 2019). By either incorporating a 45 permafrost model into the ecosystem model (Zhuang et al., 2001;Wania et al., 2009;Lyu and Zhuang, 2018) or equipping the soil model with vegetation dynamics and carbon processes , the temporal dynamics of soil temperature, permafrost dynamics and vegetation and carbon dynamics can be simultaneously simulated over cold region ecosystems. Moreover, the incorporation of detailed vadose zone and land surface processes (e.g., soil hydrology and snow cover) usually improves the model performance 50 (Lyu and Zhuang, 2018) and facilitates model ability to investigate the ecosystem response to variations in climatic and environmental conditions at various spatial-temporal scales . The importance of non-growing-season processes (e.g., freeze-thaw cycle, snow cover) was highlighted when interpreting the carbon budget observations and can significantly alter the carbon cycling and future projection of cold region ecosystems (Zhuang et al., 2001;Wania et al., 2009;Lyu and Zhuang, 2018;Zhang et al., 2018).

55
However, in most of the current modelling research in cold region ecosystems, the water and heat transfer process in the vadose zone remains independent and not fully coupled. Such consideration of vadose zone physics might result in unrealistic physical interpretations, especially for soil freezing/thawing processes (Hansson et al., 2004). In this regard, researchers have stressed the necessity to simultaneously couple the water and heat transfer process in dry/cold seasons (Scanlon and Milly, 1994;Bittelli et al., 2008;Zeng et 60 al., 2009a;Zeng et al., 2009b;Yu et al., 2016;Yu et al., 2018). Concurrently, researchers developed dedicated models, e.g., SHAW (Flerchinger and Saxton, 1989), HYDRUS (Hansson et al., 2004), MarsFlo (Painter, 2011) and its successor Advanced Terrestrial Simulator (Painter et al., 2016), and STEMMUS-FT (Yu et al., 2018;Yu et al., 2020), implementing the soil water and heat coupling physics for frozen soils (see for reviews of the relevant models in Kurylyk and Watanabe, 2013;Grenier et al., 2018;Lamontagne-Halle et al., 2020).
Promising simulation results have been reported for the soil hydrothermal regimes. While these efforts mainly focus on understanding the surface and subsurface soil water and heat transfer process (Yu et al., 2018;Yu et al., 2020) and stress the role of physical representation of freezing/thawing process (Boone et al., 2000;dynamics. 70 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 (Liu and Chen, 2000;Cheng and Wu, 2007;Yao et al., 2019). Monitoring and projecting the dynamics of hydrothermal and ecohydrological states and their responses to climate change on the Tibetan Plateau is important to help shed light on future ecosystem responses in this region. Considerable land-surface and vegetation changes have been reported in this region, 75 e.g., degradation of permafrost and variations in seasonally frozen ground thickness (Cheng and Wu, 2007;Yao et al., 2019), advancing vegetation leaf onset dates (Zhang et al., 2013), and enhanced vegetation activity at the start of growing season (Qin et al., 2016). However, there are divergences with regard to the expected ecosystem changes across the Tibetan plateau (Cheng and Wu, 2007;Zhao et al., 2010;Qin et al., 2016;. In response to climate warming, the degradation of frozen ground can positively affect 80 the vegetation growth in mountainous region (Qin et al., 2016), but it can also lead to degradation of grasslands (Cheng and Wu, 2007), depending on soil hydrothermal regimes and climate conditions (Qin et al., 2016;Wang et al., 2016).
In this study, we investigated the consequences of considering coupled water and heat transfer processes on land-surface fluxes and ecosystem dynamics in the extreme environmental conditions of the Tibetan plateau, 85 relying on land-surface and ecohydrological models confronted with multiple field observations. The inclusion or exclusion of different soil physical processes, i.e., explicitly considering the effect of soil ice content on hydrothermal properties and the tightly coupled water and heat transfer, in such environment frames the scope here. Specifically, the leading questions of the research are: i) How do different representations of frozen soil and coupled water and heat physics affect the simulated ecohydrological 90 dynamics of a Tibetan plateau meadow? ii) How does different vadose zone physics affect our interpretation of mass, energy, and carbon fluxes in the ecosystem? Answering these two questions enables evaluation of the adequacy of models in simulating feedbacks among processes and ecosystem changes across the Tibetan plateau.
In order to achieve the aforementioned goals, the detailed soil mass and energy transfer scheme developed 95 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 were fully coupled to further facilitate the model's capability in dealing with complex vadose zone processes.

Experimental Site
In this study, we make use of the Maqu soil moisture and soil temperature (SMST) monitoring network (Su et al., 2011;Dente et al., 2012;Su et al., 2013;Zeng et al., 2016), which 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. The 105 climate can be characterized by wet rainy summers and cold dry winters. The mean annual air temperature 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 a 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 about 18% organic matter for the upper soil layers (Dente et al., 110 2012;Zheng et al., 2015a;Zheng et al., 2015b;Zhao et al., 2018). The groundwater level of the grassland area fluctuates from about 8.5 m to 12.0 m below the ground surface.
For the Maqu SMST monitoring network, SMST profiles are automatically measured by 5TM 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 Planetary Boundary 115 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. Instrumentations for measuring four-component down and upwelling solar and thermal radiation (NR01-L, Campbell Scientific, Inc., USA), and liquid precipitation (T200B, Geonor, Inc., USA) are also deployed.
For this research, data from March 2016 to August 2018 collected at the central experimental site 120 (33°54'59"N, 102°09'32", elevation: 3430m) were utilized (see Figure 1). Seasonally frozen ground is characteristic of this site, with the maximum freezing depth approaching around 0.8 m under current climate conditions. The dedicated SMST profile (central station, Figure 1), with sensors installed at depths of 2.5 cm, 5 cm, 10 cm, 20 cm, 40 cm, 60 cm, and 100 cm, was used for validating the model simulations. Note that there are data gaps (25 th Mar -8 th June, 2016; 29 th Mar -27 th July, 2017, extended to 12 th Aug, 2018 for 40 125 cm) due to the malfunction of instruments and the difficulty to maintain the network under such harsh environmental conditions.

Land Surface Energy and Carbon Fluxes and Vegetation Dynamics
Starting from the raw NEE (Net Ecosystem Exchange) and ancillary meteorological data (friction velocity data (Reichstein et al., 2005). Three cases were identified according to the availability of , , and : Furthermore, we downloaded MCD15A3H (Myneni et al., 2015) and MOD17A2H (Running et al., 2015) products for this site as the auxiliary ecosystem carbon and vegetation dynamics data, from the Oak Ridge

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

T&C Model (unCPLD)
The Tethys-Chloris model (T&C) (Fatichi et al., 2012b) simulates the dynamics of energy, water, and vegetation and has been successfully applied to a very large spectrum of ecosystems and environmental conditions (Fatichi and Ivanov, 2014;Fatichi et al., 2016b;Pappas et al., 2016;Fatichi and Pappas, 2017;Mastrotheodoros et al., 2017). The model simulates the energy, water, and carbon exchanges between the 175 land surface and the atmospheric surface layer accounting for aerodynamic, undercanopy, and leaf boundary layer resistances, as well as for stomatal and soil resistance. The model further describes vegetation physiological processes including photosynthesis, phenology, carbon allocation, and tissue turnover.
Dynamics of water content in the soil profile in the plot-scale version are solved using the one-dimensional (1-D) Richards equation. Heat transfer in the soil is solved by means of the heat diffusion equation. Soil heat 180 and water dynamics are uncoupled (however, note that T&C is termed unCPLD to distinguish it later with the coupling with STEMMUS). The detailed model description is provided in the above-mentioned references and some key elements applied for this study are explained in the following.
T&C model uses the 1-D Richards equation, which describes the water flow under gravity and capillary forces in isothermal conditions 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 coordinate; S (kg m −3 s -1 ) is the sink term for transpiration and evaporation 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 190 equation to a system of ordinary differential equations in time, which can be expressed as: where , (m) is the thickness of layer i; (m s -1 ) is the vertical outflow from a layer i; Tv (m s -1 ) is the transpiration fluxes from the vegetation; , is the fraction of root biomass contained in soil layer i; Ebare (m s -1 ), evaporation from the bare soil; Es (m s -1 ), evaporation from soil under the canopy.
The heat conservation equation used in the T&C neglects the coupling of water and heat transfer physics and 195 only the heat conduction component is considered, which can be expressed as: where ρsoil (kg m −3 ) is the bulk soil density; Csoil (J kg −1 K −1 ) is the specific heat capacities of bulk soil; λeff (W m −1 K −1 ) is the effective thermal conductivity of the soil. T (K) is the soil temperature. When soil undergoes freezing/thawing process, the latent heat flux due to water phase change becomes important, which is not considered in the original T&C model, but it is in the T&C-FT (freezing/thawing) model.

T&C-FT Model (unCPLD-FT)
To account for frozen soil physics, T&C-FT model considers ice effect on hydraulic conductivity, thermal conductivity, heat capacity, and subsurface latent heat flux. However, the vapor flow and the thermal effect on water viscosity are not considered in T&C-FT, and during the non-frozen period, soil water and heat are still independently transferred as in T&C (this version is named here unCPLD-FT). To explicitly account for 205 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 θice is a prognostic variable, which is simulated along with liquid water content for each soil layer.
Specifically, when Eq. (4) is rewritten in terms of an apparent volumetric heat capacity Capp (Hansson et al., 2004;Gouttevin et al., 2012), it can be solved equivalently to Eq. (3): where Capp can be computed knowing the temperature T (K), latent heat of fusion Lf and the differential (specific) water capacity dθ/dψ at a given liquid water content θ (Hansson et al., 2004): The effective thermal conductivity λeff (W m −1 K −1 ) and the specific soil heat capacity Csoil (J kg −1 K −1 ) are computed accounting for solid particles, water, and ice content (Johansen, 1975;Farouki, 1981;Lawrence et al., 2018;Yu et al., 2018). The soil freezing characteristic curve providing the liquid water potential in a 215 frozen soil is computed following the energy conservative solution proposed by Dall'Amico et al. (2011) and it can be combined with various soil hydraulic parameterizations including van Genuchten and Saxton and Rawls to compute the maximum liquid water content at a given temperature and consequently ice and liquid content profiles at any time step (Fuchs et al., 1978;Yu et al., 2018).
Finally, saturated hydraulic conductivity is corrected in the presence of ice content (e.g., Hansson et al., 2004; 220 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 not considered here, and heat and water fluxes are still not entirely coupled in T&C-FT.

STEMMUS Model
Simultaneous Transfer of Energy, Mass and Momentum in Unsaturated Soil (STEMMUS) model solves soil 225 water and soil heat balance equations simultaneously in one time step (Zeng et al., 2011a, b;Zeng and Su, 2013). The Richards' equation with modifications made by Milly (1982) is utilized to mimic the coupled soil mass and energy transfer process. The vapor diffusion, advection, and dispersion are all taken into account as water vapor transport mechanisms. The root water uptake process is regarded as the sink term of soil water and heat balance equations, building up the linkage between soil and atmosphere (Yu et al., 2016). In 230 STEMMUS, temporal dynamics of three phases of water (liquid, vapor and ice) are explicitly presented and simultaneously solved by spatially discretizing the corresponding governing equations of liquid water flow and vapor flow.
where ρV and ρice (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 235 liquid water flow driven by the gradient of soil matric potential and temperature , respectively. ℎ and (kg m -2 s -1 ) are the soil water vapor fluxes driven by the gradient of soil matric potential and temperature , respectively. DTD (kg m -1 s -1 K -1 ) is the transport coefficient of the adsorbed liquid flow due to temperature gradient; DVh (kg m -2 s -1 ) is the isothermal vapor conductivity; and DVT (kg m -1 s -1 K -1 ) is the thermal vapor diffusion coefficient.

240
STEMMUS takes into account different heat transfer mechanisms, including heat conduction ( ), convective heat transferred by liquid and vapor flow, the latent heat of vaporization ( 0 ), the latent 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 and 245 Cice (J kg −1 K −1 ) are the specific heat capacities of soil solids, liquid, water vapor and ice, respectively; Tref (K) is the arbitrary reference temperature; L0 (J kg −1 ) is the latent heat of vaporization of 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 and qV (kg m -2 s -1 ) are the liquid and vapor water flux, respectively. Additional details 250 on the equations for solving the coupled water and heat equations can be found in Zeng et al. (2011a, b) and Zeng and Su (2013). In the appendix, a notation table is summarized for the above equations.

Coupling T&C and STEMMUS (CPLD)
As mentioned above (section 3.1-3.2), T&C considers soil water and heat dynamics independently, and T&C-FT only considers ice effects associated with latent heat, thermal and hydraulic parameters, while all other 9 soil physics processes of STEMMUS are not considered. On the other hand, while STEMMUS model can well reproduce the soil water and heat transfer process in frozen soil, it lacks a detailed description of landsurface processes and of the ecohydrological feedback mechanisms. To take advantage of the strengths of both models, we coupled STEMMUS model with the land-surface and vegetation components of T&C model (termed as CPLD) to better describe the soil-plant-atmosphere continuum (SPAC) in cold regions.

260
The current coupling procedure between STEMMUS and the T&C model is based on a sequential coupling via exchanging mutual information within one time step (see Figure 3).

Numerical Experiments
To investigate the role of increasing complexity of vadose zone physics in ecosystem functioning, three numerical experiments were designed on the basis of the aforementioned modeling framework (Table 1).
First experiment, the T&C original model was run as stand-alone, termed as unCPLD simulation. For the 280 unCPLD model, soil water and heat transfer is independent with no explicit consideration of soil ice effect.
The second experiment, the updated T&C model with explicit consideration of freezing/thawing process was run as it can estimate the dynamics of soil ice content and the related effect on water and heat transfer (e.g., blocking effect on water flow, heat release/gain due to phase change) but otherwise being exactly equal to T&C original model. This second simulation is named the unCPLD-FT simulation, where the term unCPLD 285 generally refers to the fact that T&C model and STEMMUS model are not yet coupled. The third experiment, STEMMUS model was coupled with T&C model to enable not only frozen soil physics but also additional processes and most importantly the tight coupling of water and heat effects. This simulation is named CPLD simulation. In this third scenario, vapor flow, which links the soil water and heat flow, is explicitly considered.
In addition to the ice blocking effect as presented in unCPLD-FT, the thermal effect on water flow is also 290 expressed with the temperature dependence of hydraulic conductivity and matric potential. Furthermore, not only the latent heat due to phase change, but also the convective heat due to liquid/vapor flow is simulated.
Hourly meteorological forcing (including downwelling solar and thermal radiation, precipitation, air temperature, relative humidity, wind speed, atmospheric pressure) was utilized to drive the models. For the adaptive time step of STEMMUS simulation, the linear interpolation between two adjacent hourly 295 meteorological measurements was used to generate the required values at every second. The hydrological related initial states, e.g., initial snow water equivalent, soil water and temperature profiles, were taken as close as possible to the observed ones. Since the current initial conditions of the carbon and nutrients pools in the soil are unknown, we spin-up carbon and nutrient pools running only the soil-biogeochemistry module for 1000 years using average climatic conditions and prescribed litter inputs taken from preliminary 300 simulations. Then we used the spun-up pools as initial conditions for the hourly-scale simulation over the period for which hourly observations are available. This last operation is repeated two times, which allows reaching a dynamic equilibrium of nutrient and carbon pools in the soil and vegetation.
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 305 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 literature and expert knowledge (see a summary of the adopted vegetation parameters in the supplement Table S2). All three numerical experiments shared the same soil and vegetation parameter settings.

Soil Moisture and Soil Temperature Simulations
The capability of the three models to reproduce the temporal dynamics of soil moisture is illustrated in Figure   6. By explicitly considering soil ice content, the unCPLD-FT and CPLD model captured well the response of soil moisture dynamics to the freeze-thaw cycles, while the unCPLD model lacked such capability and maintained a higher soil water content throughout the winter period, but slightly lower water content in the 340 growing season. For all three models, the consistency between the measured and simulated soil water content at five soil layers was satisfactory during the growing season, indicating the models' capability in portraying the effect of precipitation and root water uptake on the soil moisture conditions.

Soil Ice Content and Water Flux
The time-series of soil ice content and water flux from unCPLD, unCPLD-FT and CPLD model simulations for soil layers below 2 cm are presented in Figure 8. As soil ice content measurements were not available, the freezing front propagation inferred from the soil temperature measurements was employed to qualitatively assess the model performance. The phenomenon that a certain amount of liquid water flux moves upwards 365 along with the freezing front can be clearly noticed for both the unCPLD-FT and CPLD model simulations.
As the soil matric potential changes sharply during the water phase change, 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;Yu et al., 2020). Cryosuction is much more accentuated 370 in the unCPLD-FT simulation, while it is of course absent in the unCPLD model simulations (Figure 8c).
Precipitation induced downward water flux can be observed in all models during summer with very similar patterns. It is to note that compared to unCPLD-FT model, CPLD model presented a relatively lower presence of soil ice content, while its temporal dynamics was closer to the observed freezing/thawing front propagation.
The difference between the two simulations can be attributed to the constraints imposed by the 375 interdependence of liquid, ice, and vapor in the soil pores that is considered only in CPLD model.

Simulations of Land Surface Carbon Fluxes
The eddy covariance derived vegetation productivity and remote sensing ( Taking eddy-covariance observations as the reference, the onset date of grassland appears to be well captured 385 by both unCPLD and CPLD model simulations, while a delayed onset date by unCPLD-FT model. Leaf senescence and dormancy phase are a bit delayed in the models when compared with eddy-covariance data and considerably delayed when compared to MODIS-LAI, even though the latter is particularly uncertain as described above. Although there is an observable underestimation of GPP compared to the eddy covariance measurements, the dynamics of GPP, which is mainly constrained by the photosynthetic activity and 390 environmental stresses, is reasonably reproduced by all model simulations. The underestimation of GPP has magnified consequences in terms of reproducing NEE dynamics by unCPLD and CPLD models. While this might be seen as a model shortcoming, there are a number of reasons that lead to questioning the reliability of the magnitude of carbon fluxes measurements at this site. By checking other ecosystem productivity under similar conditions, the annual average GPP for the Tibetan 13 plateau meadow ecosystem ranges from 300 to 935 g C m -2 yr -1 , while the annual average NEE ranges from -79 to -213 g C m -2 yr -1 (see the literature summary in the Supplement Table S3). 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 existing estimates of the carbon exchange for such ecosystem type and are unlikely to be correct in absolute magnitude. The ecosystem respiration ( ),

400
indicating the respiration of activity of all living organisms in an ecosystem is shown in Figure 9d. 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, which further suggests the discrepancy in observed/simulated GPP is the driver of the disagreement in NEE.
The difference in the soil liquid water and temperature profile simulations between the CPLD and unCPLD lower GPP in the unCPLD simulations is instead related to a slightly enhanced water-stress induced by the different soil-moisture dynamics during the winter and summer season with a lower root zone moisture produced by the unCPLD model ( Figure 6), which affects the plant photosynthesis and growth. Differences in soil temperature profiles can also affect root respiration in generating additional small differences in GPP.

415
The energy balance closure problem, usually identified because the sum of latent (LE) and sensible (H) heat fluxes is less than the available energy (Rn-G0), is quite common in eddy covariance measurements (Su, 2002;Wilson et al., 2002;Leuning et al., 2012). The energy imbalance of EC measurements is particularly significant at sites over the Tibetan Plateau (Tanaka et al., 2003;Yang et al., 2004;Chen et al., 2013;Zheng et al., 2014). Figure 10 presents the energy imbalance of hourly LE and H by the eddy covariance 420 measurements, observed Rn by the four-component radiation measurements, and the estimated ground heat flux (G0) by CPLD model. The sum of measured LE and H was significantly less than Rn, with the slope of LE+H versus Rn equal to 0.59 (Figure 10a). Usually, the measurements of radiation are reliable (Yang et al., 2004). If we assume that the turbulence fluxes (LE, H) measurements are accurate, then the rest of energy (around 41% of Rn) should be theoretically consumed by ground heat flux G0, which is clearly impossible.

425
When compared to the available energy (Rn-G0), the slope was increased to 0.70 (Figure 10b). Table 2 demonstrated that the energy imbalance problem was 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.
14 These problems clearly suggest that care should be taken to the data mutual corroboration issue. Nevertheless, such issue is not affecting the comparison results among models with different vadose zone physics, since we did not force any parameter calibration or data-fitting procedure, but simply rely on physical constraints, literature, and expert knowledge to assign model parameters.

435
The effect of different model versions on soil water budget components is illustrated in Figure 11. The effect of coupled water and heat physics (unCPLD vs. CPLD model) on the water budget components can be summarized as: i) the amount of ecosystem water consumption ET was reduced, due to the damped surface evaporation process (evaporation from the soil surface and intercepted water). ii) water storage 460 amount 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 (see Table 1 and Supplement S1). Specifically, the high accumulation of ice content in the unCPLD-FT simulations indicates a relatively stronger cryosuction effect than in CPLD simulations. This cryosuction effect is mitigated in the fully coupled model because of water vapor transfer and thermal 465 gradients, even though different solutions in the parameterization of bulk soil thermal conductivity and volumetric soil heat capacity could also be responsible for the difference (Yu et al., 2018;Yu et al., 2020).
Overall, taking into account the fully coupled water and heat physics modifies 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). water fluxes, it is during the frozen period that water and heat transfer processes are tightly coupled (Hansson et al., 2004;Yu et al., 2018;Yu et al., 2020). Both the explicit consideration of soil ice and coupled water 490 and heat physics can affect the vadose zone water flow via altering the hydraulic conductivity and soil water potential gradients. This is testified by the fact that 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 dynamics of ice content (Figure 8), which impacts leaf onset and to a less extent hydrological fluxes. However, in the rest of the year, the simplified solution of vadose zone physics of T&C leads to very similar results as 495 the coupled one, suggesting that most of the additional physics does not modify substantially the ecohydrological response during unfrozen periods.

Conclusion
The detailed vadose zone process model STEMMUS and the ecohydrological model T&C were coupled to investigate the effect of various model representations in simulating water and energy transfer and seasonal 500 ecohydrological dynamics over a typical Tibetan meadow. The results indicate that the original T&C model tended to overestimate the variability and magnitude of soil temperature during the freezing period and the freezing-thawing transition period. Such mismatches were ameliorated by the inclusion of soil ice content and freezing-thawing processes to the original model, and further improved with explicit consideration of coupled water and heat physics. For the largest part of the simulated period (i.e., unfrozen), we found that a 505 simplified treatment of vadose zone dynamics is sufficient to reproduce satisfactory energy, water and carbon fluxes -given the uncertainty in the eddy-covariance observations. Additional complexity in vadose zone representation is mostly significant during the freezing and thawing periods as ice content simulations differ among models and the amount of water moving towards the freezing front was differently simulated. These discrepancies have an impact (even though limited to the beginning of the growing season) on vegetation 510 dynamics. The leaf onset is better captured by the unCPLD and CPLD models, while a delayed onset date was reproduced by unCPLD-FT model. Nonetheless, overall patterns for the rest of the year do not differ considerably among simulations, which suggests that the difference in vadose zone dynamics, by using a fully coupled water-heat model treatment, is not enough to affect the overall ecosystem response. This also suggests that the additional complexity might be more needed for specific vadose zone studies and 515 investigation of permafrost thawing rather than for ecohydrological applications. Nevertheless, the coupled model can reveal the hidden physically-based processes and mechanisms in the vadose zone that cannot be explained by uncoupled models, which can assist the comprehensive physical interpretations of ecosystem responses to subtle climatic changes/trends over high-altitude cold regions. In summary, our investigations using different models of vadose zone physics can be helpful to support the development and application of 520 earth system models as they suggest that a certain degree of complexity might be necessary for specific analyses.   Independent water and heat transfer: Soil water and heat transfer process is independent.
FT induced water and heat transfer coupling: Soil water and heat transfer process is coupled only during the freezing/thawing (FT) period. Soil water flow is affected by temperature only through the presence of soil ice content 815 (the impedance effect).
Tightly coupled water and heat transfer: Soil water and heat transfer process is tightly coupled; vapor flow, which links the soil water and heat flow, is taken into account; thermal effect on water flow is considered (the hydraulic conductivity and matric potential is dependent on soil temperature; when soil freezes, the hydraulic conductivity is reduced by the presence of soil ice, which is temperature dependent); the convective/advective heat due to liquid/vapor 820 flow can be calculated.
Ice effect on soil properties: the explicit simulation of ice content and its effect on the hydraulic/thermal properties.