Sensitivity Analysis of the Noah-MP Land Surface Model for Soil Hydrothermal Simulations Over the Tibetan Plateau

Abstract

The Tibetan Plateau (TP) features unique and highly heterogeneous soils, terrains, vegetation, and climate. Accurately modeling complex freeze-thaw processes and their hydrothermal impacts remains a great challenge. This study focused on deciphering the spatiotemporal variability of diverse parameterization schemes in the soil hydrothermal simulations using the Noah-MP land surface model. We first discussed the spin-up time required by the model to reach the equilibrium state, and then performed a sensitivity analysis of these schemes. The Moderate Resolution Imaging Spectroradiometer land surface temperature and Soil Moisture Active Passive remote sensing products were used as benchmarks to evaluate the schemes' performance. Results show that longer spin-up times are required in permafrost regions owing to water phase changes. Ground temperature and soil temperature are mainly sensitive to energy-related schemes. Vegetation-related schemes play an important role after the growing season begins on the southeastern TP. Soil water content shows strong sensitivity to schemes related to both water and energy transport. However, the sensitivity of these energy-related schemes is weakened when simulating total soil moisture, including the total amount of water and ice, indicating that these schemes have marked impacts on soil freeze-thaw processes. These results reveal the different spatial (both regional and depth-related) and temporal effects of parameterization schemes; we also provided a preliminary selection of these schemes at a regional scale that could facilitate the further improvement of the soil hydrothermal simulations on the TP.

Introduction

Land surface processes play an important role in the climate system by controlling land-atmosphere exchanges of energy, momentum, and mass (Dickinson et al., 2006; Hurrell et al., 2013; Pitman, 2003). Many numerical and physics-based land surface models (LSMs) have been developed to understand the physical mechanisms of land surface processes (F. Chen & Dudhia, 2001; Dai et al., 2003; Niu et al., 2011; Oleson et al., 2008; Sellers et al., 1986), and are further used in weather and climate prediction (Dickinson et al., 1986; Koster et al., 2006; Sellers et al., 1996). However, there are still considerable uncertainties and deficiencies within land surface simulations (Dirmeyer et al., 2006; Gayler et al., 2013; Koster et al., 2006), which mainly result from input data, model parameters, and model structure (Clark et al., 2011; Pastres & Ciavatta, 2005; Radwan et al., 2004; G. Zhang et al., 2016). Model structure uncertainty is commonly caused by spatiotemporal discretization or parameterization schemes of physical processes (Parrish et al., 2012). Several studies used multimodel and multischeme methods to investigate the structural uncertainty of LSMs (Gayler et al., 2014; Z.-L. Yang et al., 2011; Zhou et al., 2012), and dedicated projects were launched using such methods; these include the Project for Intercomparison of Land Surface Parameterization Schemes (PILPS, Henderson-Sellers et al., 1993) and the Global Soil Wetness Project (GSWP, Dirmeyer, 2011). Such research has shown that disparities in parameterization schemes cause differences in LSM behaviors; that is, even when driven by the same meteorological forcing, large discrepancies remain between various LSM simulations (Jimenez et al., 2011; Xia et al., 2014; H. Zheng et al., 2019); this has hampered the understanding of complex land surface processes. Therefore, model simulations require comprehensive evaluation, and the ways in which disparities are caused within parameterization schemes need to be further investigated.

Individual influence is always buried in the integration of the effects, which means that previous efforts have not fulfilled the objective of systematically quantifying uncertainties in parameterization schemes (Clark et al., 2011). Through sensitivity analysis of model outputs, which is an effective method to identify informative parameterizations and select model parameters and schemes, the physical mechanisms of land surface processes may be discovered and the simulation performance of LSMs may be improved. In previous work, the sensitivity of different parameters (L. Chen et al., 2016; Cuntz et al., 2015; Franks & Beven, 1997; Pitman, 1994) and different schemes (Cai et al., 2014; Gan et al., 2019; J. Li et al., 2022; X. Li et al., 2020; H. Zheng et al., 2019) has been comprehensively assessed, and it has been demonstrated that the sensitivity varies with the study location, model structure, input data, and methods used (Razavi & Gupta, 2015; Rosero et al., 2010). In reality, a combination of schemes that performs well at one site may not be appropriate for another site, owing to heterogeneity in site characteristics and their associated physical processes (Clark et al., 2011; Gan et al., 2019; Gayler et al., 2014; Mendoza et al., 2015). Without proper constraints, existing models can show a large range in performance level, and improper selection of parameterization, together with inaccuracies in model formulation, will likely produce substantial model errors (Godfrey & Stensrud, 2010; van den Hurk et al., 2011). Therefore, by combining multiple sites and remote sensing data, regional assessments of sensitivity are required to improve model performance and discover the underlying relationships between model outputs and other characteristics (Gan et al., 2019; N. Ma et al., 2017; Wang et al., 2018; H. Zheng & Yang, 2016).

Soil hydrothermal processes related to energy and water transfer are key factors in the control of land surface simulations. In particular, soil freeze-thaw processes can bring about simulation differences by changing soil water and heat transfer (Nicolsky et al., 2007; Niu & Yang, 2006; Qin et al., 2016; Vecellio et al., 2019; Walvoord & Kurylyk, 2016; Woo et al., 2008). The Tibetan Plateau (TP) or “Third Pole” (Qiu, 2008) is characterized by complex terrain and heterogeneous land cover and soil structure, and has a widely distributed frozen soil; the total area of permafrost, seasonally frozen ground, and unfrozen ground on the TP is 1.06 × 10⁶ (40%), 1.46 × 10⁶ (56%), and 0.03 × 10⁶ km² (1%), respectively (Zou et al., 2017). Studies have shown that the simulated energy and water fluxes over the TP are dominated by soil parameterization schemes (Y. Chen et al., 2012; Gao et al., 2015; X. Li et al., 2020; Lu et al., 2020; Luo et al., 2017; Song et al., 2020; D. Zheng et al., 2017). However, many of the LSMs have difficulties, for example, in reliably representing the soil temperature (ST) and top-layer soil moisture, which have been systematically underestimated in the TP (Bi et al., 2016; K. Yang et al., 2009), and marked discrepancies have been observed in the responses of the LSMs to changes in soil parameterization (J. Li et al., 2018; Q. Li et al., 2021; X. Li et al., 2021; Pan et al., 2017). Recent efforts have improved the representation of soil physics in LSM simulations in this region, which has been achieved in three respects. First is the implementation of improved parameterization schemes such as improving ST and moisture profiles by implementing new under-canopy turbulence and root water uptake strategies (D. Zheng et al., 2015) and a scheme considering diurnal variation of thermal roughness length (Y. Chen et al., 2010; K. Yang et al., 2008). Second is consideration of the heterogeneity of the soil profile, including stratified active layers (Pan et al., 2017), soil organic matter, and gravel (Y. Chen et al., 2012; Y. Liu et al., 2021; Y. Sun et al., 2017; Yuan et al., 2021). Third is the introduction of new physical processes such as rhizosphere physics in the topsoil (Gao et al., 2015), ground surface deformation (S. Liu et al., 2022), and the modification of the distinction between freezing and thawing by using the virtual temperature instead of a constant freezing point (Q. Li et al., 2021). Nonetheless, owing to the harsh environment of the TP, there is insufficient land surveying and a lack of observation stations to provide data for the sensitivity analysis of LSMs across the entire region (Shi et al., 2004; Suh & Lee, 2004). Hence, most previous studies have been site-specific and their applicability to the whole region remains to be tested. Fortunately, advanced satellite observations from space have greatly facilitated LSMs evaluations and calibrations at regional or global scales (Anav et al., 2015; Toure et al., 2016; Xia et al., 2012).

The main purpose of this paper is to discover the source of uncertainty in soil hydrothermal simulations (in different seasons and at different depths) across the entire TP based on multi-parameterization schemes and then to use remote sensing data as “ground truth” data to select suitable options for these schemes. In Section 2, we describe the Noah-MP LSM, data, experimental design, and methods used in the study. We present our results in Section 3 and our conclusions in Section 4.

Descriptions of Model, Data, Experimental Design, and Methods Noah-MP Land Surface Model

We used Noah-MP LSM version 4.1 in this study. The Noah-MP was developed from the Noah LSM and enhanced the description of several physical processes. It provides a unique multi-parameterization framework that allows the model to be run with different combinations of alternative schemes for individual land processes (Niu et al., 2011), and it can facilitate the study of large-scale parameterization sensitivity. The user can choose from a variety of schemes to improve the simulations and identify the physical mechanisms underlying the results (Cai et al., 2014; Gan et al., 2019; Hong et al., 2014; Hu et al., 2020; X. Li et al., 2021; Niu et al., 2011; Z.-L. Yang et al., 2011; You et al., 2020; H. Zheng et al., 2019), and further combine this model with Weather Research and Forecast Model for weather forecast and short-term climate predictions (Barlage et al., 2015). These schemes were listed in Table 1; each process contains 2–4 different options to select, which correspond to different parameterization schemes (We briefly compared the different options in Text S1 in Supporting Information S1).

Table 1 Noah-MP Parameterization Schemes and Their Options Used in This Study

Parameterization schemes	Options
DVEG: Dynamic vegetation	1. Off (preset leaf area index and vegetation fraction)
	2. On (together with CRS = 1)a
CRS: Canopy stomatal resistance	1. Ball-Berry: Relates stomatal conductance to CO2 exchange rate through stomata (Ball et al., 1987)a
	2. Jarvis: As a function of soil moisture, atmospheric temperature, radiation availability, and vapor pressure deficit (Jarvis, 1976)
BTR: Soil moisture factor for stomatal resistance	1. Noah: Function of soil moisture, as in Noah (F. Chen & Dudhia, 2001)a
	2. CLM: Matric potential related, as in CLM (Oleson et al., 2004)
	3. SSiB: Matric potential related, as in SSiB (Y. Xue et al., 1991)
RUN: Runoff and groundwater	1. TOPMODEL-based runoff with simple groundwater (Niu et al., 2007)a
	2. TOPMODEL-based runoff with an equilibrium water table (Niu et al., 2005)
	3. Noah infiltration-excess-based surface runoff with free drainage (Schaake et al., 1996)
	4. BATS runoff with free drainage (Z.-L. Yang & Dickinson, 1996)
SFC: Surface layer drag coefficient	1. M-O: Monin-Obukhov similarity theory (Brutsaert, 1982)a
	2. Chen97: Neglecting the zero-displacement height (F. Chen et al., 1997)
FRZ: Supercooled liquid water	1. NY06: No iteration, generalized freezing point depression (Niu & Yang, 2006)a
	2. Koren99: Koren's iteration, variant freezing point depression (Koren et al., 1999)
INF: Frozen soil permeability	1. NY06: Linear effects, more permeable, function of soil moisture (Niu & Yang, 2006)a
	2. Koren99: Nonlinear effect, less permeable, function of liquid water volume (Koren et al., 1999)
RAD: Radiation transfer	1. Modified two-stream, canopy gaps from 3D structure and solar zenith angle (Niu & Yang, 2004)a
	2. Two-stream applied to grid-cell, no canopy gap (Niu & Yang, 2004)
	3. Two-stream applied to vegetated fraction, gaps from vegetated fraction (Niu & Yang, 2004)
ALB: Ground snow surface albedo	1. BATS: Computes snow surface albedo for direct and diffuse radiation over visible and near-infrared wave bands (Z.-L. Yang et al., 1997)
	2. CLASS: Computes the overall snow surface albedo accounting for fresh snow albedo and snow age (Verseghy, 1991)a
SNF: Partitioning precipitation into rainfall or snowfall	1. Jordan: Relatively complex functional form (Jordan, 1991)a
	2. BATS: Assumes all precipitation as snowfall when air temperature below freezing temperature plus 2°
	3. Assumes all precipitation as snowfall when air temperature below freezing temperature
TBOT: Lower boundary condition of soil temperature	1. Zero heat flux from bottom
	2. TBOT at 8 m depth read from a preset filea
STC: Snow and soil temperature time scheme	1. Semi-implicit, flux top boundary conditiona
	2. Noah: Fully implicit, temperature top boundary condition

Note. CLM is Community Land Model; SSiB is Simplified Simple Biosphere Model; BATS is Biosphere-Atmosphere Transfer Scheme; CLASS is Canadian Land Surface Scheme; and TOPMODEL is Topography-based hydrological model.

^aDefault option.

Data Used in the Study Model Input

The forcing data determines the accuracy of simulations to a certain extent in off-line land surface modeling. For example, Lu et al. (2020) validated six forcing data sets and found that these data sets have discrepancies in representing precipitation, and the China Meteorological Forcing Data set (CMFD), developed at the Institute of TP Research, Chinese Academy of Sciences (He et al., 2020; K. Yang et al., 2010), has better consistency with observations. Therefore, this data set was used as the forcing data in this study, which includes 10-m wind speed, 2-m air temperature and humidity, surface air pressure, downward shortwave and longwave radiation, and precipitation, recorded at 3-hourly temporal resolution and $$0.1{}^{\circ}\times 0.1{}^{\circ}$$ spatial resolution, with good evaluation from previous land surface modeling studies across China (Y. Chen et al., 2011; Leng et al., 2015; B.-L. Xue et al., 2013).

This study aimed to conduct a preliminary analysis of the complex soil physical processes occurring across the TP, so it was necessary to select reliable soil texture data, as simulation conditions can have a great influence on regional simulations under different climatic regimes (H. Zheng & Yang, 2016). The default soil texture data set of Noah-MP is from the Food and Agriculture Organization. However, this data set describes the soil texture of the TP as almost entirely loam soil (Figure 1d in Gan et al. (2019)), which is not sufficiently specific to explore the complex soil physical processes that occur across this region. Therefore, we replaced this data set with a new soil property data set released by Beijing Normal University (hereafter referred to as the BNU data set, Shangguan et al., 2013). Fang et al. (2021) compared these two datasets and demonstrated that the BNU data set can better describe the soil characteristics of the TP. We used the soil loam, clay, and sand contents from this data set and soil pedo-transfer functions (Saxton & Rawls, 2006) provided in the model to convert these content values to the soil parameters required for the model simulations.

Other static input data required by the model were obtained from the Weather Research and Forecasting Model—Preprocessing System data set, and the initial field data came from the Global Land Data Assimilation System. We provided distributions of elevation, vegetation category, and permafrost in Figure 1.

Figure 1. Distributions of (a) elevation, (b) vegetation category, and (c) permafrost.

Data Used for Selecting Schemes

In this study, two satellite products were selected as “ground truth” data and compared with the simulation results of ground temperature (GT) and soil water content (SWC) in the first layer (5 cm) to select better parameterization scheme options. The corresponding moments of simulations were used because the remote sensing data were observed at a particular time; and we selected the schemes when there are more than 20% effective remote sensing values (including 20% daytime values and 20% nighttime values) in each season for analysis, which makes more grid points have scheme selections and makes the results have good reliability. For grids that don't meet this standard, we recommend using the default scheme combination for simulations.

The Moderate Resolution Imaging Spectroradiometer (MODIS) land surface temperature product was used in this study. The MOD11A1 product is produced from the Terra satellite with overpass times of approximately 10:30 a.m. (local solar time) in descending orbit and 10:30 p.m. in ascending orbit, and the MYD11A1 product is produced from the Aqua satellite with overpass times of approximately 1:30 p.m. in ascending orbit and 1:30 a.m. in descending orbit. These values were used only if they passed quality control. Previous studies have verified the accuracy of these data and have used them for model calibration and validation (Y. Chen et al., 2011; Y. Chen et al., 2017; X. Chen et al., 2017; Y. M. Ma et al., 2014; Zhong et al., 2010). The Soil Moisture Active Passive (SMAP) satellite is an L-band satellite that can provide global-scale soil moisture and freeze/thaw-state data. It was launched on 31 January 2015, and can observe the Earth's surface twice a day at sun-synchronous 6:00 a.m. in descending orbit and 6:00 p.m. in ascending orbit. Y. Chen et al. (2017) evaluated SMAP and found that it can successfully capture the amplitude and temporal variation of soil moisture at different observation stations across the TP.

Experimental Design Simulation Setup

Individual simulations in each season, rather than year-round simulations, were used to explore the influence of parameterization schemes on seasonal simulations. All pre-simulation spin-ups were performed using the default combination. We then used the restart function of the model and selected different options at the beginning of each season, which resulted in each season having separate simulations (Figure 2). This method ensured that any discrepancies were caused only by the selection of different options in a specific season and not by output data discrepancies from previous seasons. The analysis was based on hourly simulations from 2015 to 2018.

Figure 2. A flow chart of seasonal simulations.

Notably, with the development of the Noah-MP LSM, new parameterization options have been continually added to the model. To save computing resources, only the schemes provided on the Noah-MP official website (https://www.jsg.utexas.edu/noah-mp/scheme-options) were used in this study. Furthermore, the total number of possible scheme combinations stands at 4608 (Gayler et al., 2014), and as the model develops, the number reaches 2764800 in the latest version (J. Li et al., 2022); this leads to the consumption of huge amounts of computing time when considering the complex interactions and dependencies between scheme combinations over a large area using all combinations. Thus, the interactions between these schemes were not considered in this study; that is, the default options set was used for other schemes when changing an option of a specific parameterized scheme (Table 2).

Table 2 The Scheme Options Used in This Study

Experiments	Schemes
	ALB	BTR	CRS	DVEG	FRZ	INF	RAD	RUN	SFC	SNF	STC	TBOT
Default	2	1	1	2	1	1	1	1	1	1	1	2
ALB1	1	1	1	2	1	1	1	1	1	1	1	2
BTR2	2	2	1	2	1	1	1	1	1	1	1	2
BTR3	2	3	1	2	1	1	1	1	1	1	1	2
CRS2	2	1	2	1	1	1	1	1	1	1	1	2
DVEG1	2	1	1	1	1	1	1	1	1	1	1	2
FRZ2	2	1	1	2	2	1	1	1	1	1	1	2
INF2	2	1	1	2	1	2	1	1	1	1	1	2
RAD2	2	1	1	2	1	1	2	1	1	1	1	2
RAD3	2	1	1	2	1	1	3	1	1	1	1	2
RUN2	2	1	1	2	1	1	1	2	1	1	1	2
RUN3	2	1	1	2	1	1	1	3	1	1	1	2
RUN4	2	1	1	2	1	1	1	4	1	1	1	2
SFC2	2	1	1	2	1	1	1	1	2	1	1	2
SNF2	2	1	1	2	1	1	1	1	1	2	1	2
SNF3	2	1	1	2	1	1	1	1	1	3	1	2
STC2	2	1	1	2	1	1	1	1	1	1	2	2
TBOT1	2	1	1	2	1	1	1	1	1	1	1	1

Note. DVEG must be set to off when CRS is the Jarvis scheme based on the Noah-MP model's structure.

Model Spin-Up Process

To eliminate the anomalies caused by the initial input data, simulated variables should first reach an equilibrium state, and spin-up is therefore needed. The spin-up time is broadly defined as an adjustment process as the model approaches equilibrium following initial soil anomalies or after some abnormal environmental forcings. Without the spin-up process, the uncertainties can propagate from the beginning of the modeling chain to the final simulations, and the results can be misleading (Ajami et al., 2014; Cosgrove et al., 2003; Z.-L. Yang et al., 1995). Z.-L. Yang et al. (1995) indicated that land surface schemes require many years to come to thermal and hydrologic equilibrium with the forcing meteorology, and the time needed depends on the total moisture holding capacity and the initialization of the moisture stores. Cai et al. (2014) suggested that the model requires 34 years of spin-up time while simulating water table depth for the Mississippi River basin. Gao et al. (2015) pointed out that for the TP, the spin-up time varies from 4 to 30 years depending on different soil physical processes.

In this study, we performed a spin-up test to determine how long the model takes to reach equilibrium when simulating soil hydrothermal processes across the TP. Equation 1 (Cai et al., 2014; L. Chen et al., 2016; Gao et al., 2015; Z.-L. Yang et al., 1995) shows that the equilibrium state is satisfied when the variable difference between the simulations of two consecutive years becomes <0.1% of the annual mean ($$\mathrm{V}\mathrm{a}\mathrm{r}$$ is the variable, n is the years for spin-up, and the spin-up values are the ratios of the variable difference to the annual mean). However, when the spin-up time is too long in the simulations of a large area, these annual simulations will require considerable computing time. Therefore, Equation 2 is used to determine cases where the spin-up time is longer than 10 years; that is, the simulation difference between two consecutive 5-years periods become less than 0.5% of the annual mean, and thus n + 5 is the required years for spin-up. Under this equation, spin-up for more than 10 years needs to be simulated only once every 5 years. [Image Omitted. See PDF]

Statistical Methods Taylor Score

In this study, the Taylor score was used for both sensitivity analysis of parameterization schemes and the selection of suitable schemes. Taylor (2001) first proposed this statistical method, which comprehensively considers discrepancies in the correlation coefficient and standard deviation between different data in a unified interval. This method has since been further developed; for example, Hong et al. (2014) not only considered the correlation coefficient and normalized standard deviation but also added the normalized average value into the variable, creating a more comprehensive statistical variable. To combat negative correlation coefficients, Gan et al. (2019) modified the correlation coefficient term. In this study, we combined the previous modifications of the Taylor score with Equation 3, which takes into account both the differences in average value and the negative correlations. [Image Omitted. See PDF]Where, $$S$$ is the Taylor score, $$R$$ is the correlation coefficient, and $${\sigma }_{\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}}$$ and $${\upsilon }_{\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}}$$ are the ratios of the normalized standard deviation and average value of two data sets, respectively. A lower S value indicates a larger difference between the two groups of data. In Sector 3.2, the sensitivity was determined by comparing the discrepancies between the simulations using different schemes (if there were more than two options in a scheme, a pairwise comparison was performed and the average value was taken), with a lower S value indicating higher sensitivity. On the contrary, the highest S value indicates that two sets of data have the best consistency, which indicates that the simulation is closest to the remote sensing data, and this option is selected as the preferred option of the parameterization scheme in Sector 3.3.

Kolmogorov-Smirnov Test

In Section 3.2, except for the most sensitive scheme, we used a regionally averaged Taylor score across the entire TP region to rank the color bars (top to bottom correspond to sensitivity from high to low) and used the Kolmogorov-Smirnov test (K–S test) to verify the reliability of these rankings. The Kolmogorov-Smirnov test is a nonparametric hypothesis test that is mainly used to test whether a set of samples (one-sample K–S test) come from a probability distribution or to compare whether two sets of samples (two-sample K–S test) have the same distribution (Fasano & Franceschini, 1987; Justel et al., 1997; Lopes, 2011). In this study, the two-sample K–S test was applied to verify whether different options resulted in significantly different average values (red stars indicate significant differences in average values) and standard deviations (blue stars indicate significant differences in standard deviations) at the 95% confidence level for the simulations of the entire TP region.

Results Spin-Up Time for Simulating Soil Hydrothermal Processes

Figure 3 shows the time taken for GT, ST, SWC, and total soil moisture (TSM) to reach equilibrium states and pass the spin-up test at different depths across the TP.

Figure 3. Spin-up time for soil hydrothermal simulations across the Tibetan Plateau. (a–e) Are ground temperature and soil temperature at 5, 25, 70, and 150 cm, respectively; (f–i) are soil water content at 5, 25, 70, and 150 cm, respectively; and (j–m) are total soil moisture at 5, 25, 70, and 150 cm, respectively.

Figures 3a–3e, which is used to describe the temperature simulations (GT and ST), shows that different regions required different amounts of time to reach an equilibrium state from the surface to the deep layers. The northern part of the TP generally takes longer to reach equilibrium than the southern part, except for the Qaidam Basin, which suffers extreme droughts and has a relatively low elevation. This pattern is highly consistent with the spatial distribution of permafrost (Zou et al., 2017), which indicates that it takes a long time to establish equilibrium in the areas where the underlying layers comprise permafrost, and soil that is seasonally frozen finishes the spin-up in a relatively shorter time. Similar to the results from Cosgrove et al. (2003), spin-up time exhibits large spatial variability. From the vertical perspective of different depths, this phenomenon not only occurs in the deep soil where permafrost exists but also has the same distribution characteristics in the shallow layers and ground, indicating that the spin-up time of the shallow layers and ground is also affected by the equilibration of the freeze-thaw process at deep layers. At the same time, different depths require slightly different amounts of time to pass the spin-up test.

As shown in Figures 3f–3i and Figures 3j–3m, SWC reaches equilibrium slowly in the northern part of the TP where permafrost exists, while seasonally frozen soil reaches equilibrium more quickly. Except in the arid Qaidam Basin, where the model required much longer time for the water table depth to reach an equilibrium state (Niu et al., 2007), the overall time required is shorter than for temperature simulations. Meanwhile, the spin-up time required for deep layers increases compared to the shallow layers, especially in this basin. In addition, in the simulations of TSM, which considers both soil water and ice, the spin-up time is shorter than that of SWC, which only considers liquid water, indicating that the time is relatively short without considering the phase transformation between water and ice. Therefore, the freeze-thaw processes in the soil markedly increase the spin-up time in permafrost regions.

Figure 4 displays the regional average spin-up process of the entire TP. To reach equilibrium, the GT, ST, SWC, and TSM require 6, 9, 5, and 4 years, respectively. However, in this spin-up test, it takes a long time for the permafrost regions and the Qaidam Basin to reach equilibrium. As a whole, 30 years is proven to be enough time for most regions of the TP to pass the test, and this is consistent with site observations across the TP (Gao et al., 2015). Therefore, spin-up for 30 years was carried out in subsequent experiments prior to the simulations.

Figure 4. The model spin-up process for (a) ground temperature and soil temperature, (b) soil water content, and (c) total soil moisture, which pass the spin-up test when the curves fall into the green areas.

Spatiotemporal Distribution of Parameterization Sensitivity

Sensitivity analysis is an effective method to explore model mechanisms and improve their performance. To better understand the spatiotemporal distribution characteristics of the sensitivity, we analyzed the sensitivity of different soil hydrothermal variables, depths, and seasons across the TP; this allowed us to explore the physical mechanisms behind the complex variations in soil hydrothermal characteristics.

Ground Temperature and Soil Temperature

Figures 5–8 display the spatial patterns of the most sensitive schemes across the TP when simulating the GT and ST. Specific patterns emerge for different seasons, as well as for different depths and regions.

Figure 5. The most sensitive scheme of (a) ground temperature, (b) soil temperature (ST) at 5, (c) ST at 25, (d) ST at 70, and (e) ST at 150 cm over the Tibetan Plateau (TP) in winter. The rank of the color bars on the right of each panel indicates the sensitivity over the entire TP, and the sensitivity decreases from top to bottom. The red and blue stars to the right of each color bar indicate that this scheme will bring significant differences to the temperature simulations at the level of the entire TP based on the K-S test at the 95% confidence level.

Figure 6. The same as Figure 5 but in spring.

Figure 7. The same as Figure 5 but in summer.

Figure 8. The same as Figure 5 but in Autumn.

On the surface and at shallow depths (5 and 25 cm) in winter (Figures 5a–5c), the surface layer drag coefficient (SFC) and snow/soil temperature time (STC) schemes are the main components of most sensitivity schemes. The ground snow surface albedo scheme (ALB) covers the Qaidam Basin; all these schemes are highly related to energy transfer and thermal regime (Q. Li et al., 2021; Xie et al., 2017; K. Yang et al., 2008; T. Zhang, 2005). The southeastern Tibetan Plateau (TP) is mainly covered by a mixture of ALB and dynamic vegetation at this time. With increasing depth (Figures 5c–5d), the area most sensitive to ALB expands in the northeast of the TP, and vegetation-related schemes still show high sensitivity in the southeast of the TP, owning to transpiration and its cooling effect (Shen et al., 2015). As the depth increases, the main role on the western TP is taken over by the lower boundary condition of the soil temperature schemes (TBOT). The sensitivity of the supercooled liquid water scheme (FRZ) and TBOT gradually occupy the majority of the sensitivity regions in the deepest layers (150 cm, Figure 5e) on the northern and southern TPs, respectively.

In spring (Figure 6), except for the southeastern TP (including regions covered with forests and grasslands), which is most sensitive to the vegetation-related dynamic vegetation (DVEG) and canopy stomatal resistance (CRS) schemes, most areas are mainly sensitive to SFC in shallow layers and STC in deep layers.

With the growth of vegetation on the plateau in summer (Figure 7), the sensitivity of vegetation-related dynamic vegetation and CRS reaches its maximum, and the area sensitive to vegetation schemes gradually expands to the entire eastern and southern parts of the Tibetan Plateau (TP) (including regions covered with forests, grasslands, and mixed grasslands) and reaches a deeper depth than that in spring. In other areas of the northern TP (including barren or sparsely vegetated regions), the shallow layer is still sensitive to the SFC and STC schemes, while the bottom layer is converted to the TBOT scheme.

Features in autumn (Figure 8) are similar to those in spring and summer in the shallow and middle layers; however, the most sensitive areas related to vegetation schemes gradually decrease in the deep soil and are completely replaced by the TBOT in the bottom layer.

Results from the rank of color bars (Figures 5–8) and the K-S test show that the selections of schemes lead to significant differences in regionally averaged simulations of GT and ST, which mainly come from the SFC, DVEG, CRS, ALB, STC, radiation transfer (RAD), FRZ, TBOT, and runoff and groundwater (RUN) schemes in winter; the SFC, STC, DVEG, CRS, ALB, frozen soil permeability (INF), partitioning precipitation into rainfall and snowfall (SNF), RAD, FRZ, RUN, and TBOT schemes in spring; the DVEG, CRS, SFC, STC, INF, RAD, SNF, RUN, and TBOT schemes in summer; and the SFC, DVEG, CRS, ALB, STC, RAD, TBOT, FRZ, and INF in autumn. Figures 9a–9d quantifies the regionally average sensitivity in different seasons (Table S1 in Supporting Information S1 shows the specific values of the sensitivity). In the seasons with freeze-thaw processes, especially in spring, the uncertainty caused by scheme selection is more obvious. Meanwhile, sensitivity varies greatly depending on depth and season.

Figure 9. Regionally average sensitivity of ground temperature and soil temperature simulations in (a) winter, (b) spring, (c) summer, and (d) autumn. S is the Taylor Score.

Soil Water Content

Compared with GT and ST, the most sensitive schemes of SWC are more complex because of the freeze-thaw processes within the soil.

Figures 10–13 show the most sensitive schemes of SWC simulations across the TP, at different soil depths, and during different seasons. In winter, the upper soil layer (Figure 10a) shows the greatest sensitivity to vegetation-related DVEG and CRS schemes on the southeastern TP (including regions covered with forests and grasslands), freeze-thaw-related INF and FRZ schemes on the northwestern TP (including regions containing permafrost), and the runoff and groundwater scheme (RUN) on parts of the northern TP; these schemes are all directly or indirectly related to water transport (Niu et al., 2011). Other regions are mainly covered by the energy-related STC scheme. Except for the soil hydrothermal process, the STC scheme also influences the snow process, which could further affect the simulations (X. Li et al., 2021). As the soil depth increases, the supercooled liquid water scheme (FRZ) becomes the most sensitive scheme in most regions of the TP (Figures 10b–10d).

Figure 10. The most sensitive scheme of (a) soil water content (SWC) at 5, (b) SWC at 25, (c) SWC at 70, and (d) SWC at 150 cm over the Tibetan Plateau (TP) in winter. The rank of the color bars on the right of each panel indicates the sensitivity over the entire TP, and the sensitivity decreases from top to bottom. The red and blue stars to the right of each color bar indicate that this scheme will bring significant differences to the SWC simulations at the level of the entire TP based on the K-S test at the 95% confidence level.

Figure 11. The same as Figure 10 but in spring.

Figure 12. The same as Figure 10 but in summer.

Figure 13. The same as Figure 10 but in autumn.

When the frozen soil begins to melt in spring, energy-related schemes (STC and SFC) occupy most of the permafrost at shallow to middle soil depths (up to 70 cm); the exception is the southeast TP, which is sensitive to the INF scheme at shallow depths and vegetation-related schemes at middle depths (Figures 11a–11c) because SWC is affected by transpiration rates in this region with lush vegetation (Cai et al., 2014). At 150 cm depth (Figure 11d), most of the TP remains sensitive to INF because the active layer of permafrost thaws more slowly here than in the shallow and middle layers. Only parts of the southwestern TP and Qaidam Basin, located in the seasonal freezing soil region, are most sensitive to the energy-related SFC scheme at this depth. This result implies that the simulations of energy-related evaporation and freeze-thaw processes play important roles in SWC simulations across most areas of the TP in spring.

In summer (Figure 12), although the frozen soil at the shallow and middle layers melts in most parts of the permafrost region, the SWC simulations of each layer are still highly sensitive to the frozen soil-related INF scheme across the permafrost region, reflecting the response of the shallow and middle soil layers to the changes in the deep freeze-thaw processes because the two options have the same permeability at the layers without ice. In addition, the SWC simulations are sensitive to the runoff and groundwater scheme (RUN), vegetation-related schemes (DVEG and CRS), and energy-related SFC scheme, which appear in certain areas of the southeastern and western parts of the TP, the Qaidam Basin, and regions that comprise seasonal frozen and unfrozen soil.

When the soil freezes again in autumn, SWC at shallow and middle depths (Figures 13a–13c) in the central and western parts of the TP become sensitive to the energy-related SFC and TBOT schemes again, indicating that differences in the simulations of the freeze-thaw processes and evaporation result in significant discrepancies in the SWC simulations. The northern part of the TP in autumn is covered by the supercooled liquid water scheme (FRZ), while the southeastern part is still controlled by the vegetation-related DVEG scheme (including regions covered with forests and grasslands). At 150 cm depth (Figure 13d), most regions are affected by the runoff and groundwater scheme (RUN); the exception is the northwestern parts of the permafrost regions, which are most sensitive to the TBOT scheme related to the soil freezing process.

In the regionally average simulations, SWC is sensitive to STC, SFC, FRZ, DVEG, CRS, ALB, INF, RUN, RAD, and TBOT in winter; STC, INF, SFC, ALB, DVEG, CRS, FRZ, RUN, TBOT, ALB, and RAD in spring; INF, STC, SFC, RUN, CRS, DVEG, RAD, and TBOT in summer; and DVEG, CRS, SFC, STC, FRZ, INF, ALB, RUN, RAD, and TBOT in autumn. Similar to GT and ST, the sensitivity of the schemes is stronger in winter and spring, and the sensitivity of the INF scheme on all layers is greater than that of other schemes in summer (Figure 14; for the specific values of this sensitivity, see Table S2 in Supporting Information S1).

Figure 14. Regionally average sensitivity of soil water content simulations in (a) winter, (b) spring, (c) summer, and (d) autumn. S is the Taylor Score.

Total Soil Moisture Considering Water and Ice

Figures 15–18 show the sensitivity of TSM. When considering the total content of liquid water and solid ice in the soil layers, the most sensitive schemes of TSM simulations are simpler than those of SWC simulations because the part of the sensitivity caused by soil freeze-thaw processes is filtered out when the direct effects of phase change are not considered.

Figure 15. The most sensitive scheme of (a) total soil moisture (TSM) at 5, (b) TSM at 25, (c) TSM at 70, and (d) TSM at 150 cm over the Tibetan Plateau (TP) in winter. The rank of the color bars on the right of each panel indicates the sensitivity over the entire TP, and the sensitivity decreases from top to bottom. The red and blue stars to the right of each color bar indicate that this scheme will bring significant differences to the TSM simulations at the level of the entire TP based on the K-S test at the 95% confidence level.

Figure 16. The same as Figure 15 but in spring.

Figure 17. The same as Figure 15 but in summer.

Figure 18. The same as Figure 15 but in autumn.

At shallow soil depths (Figures 15a, 15b, 16a, 16b, 17a, 17b, 18a, and 18b), the scheme of frozen soil permeability (INF) occupies the highest sensitivity position throughout the year. This means that permeability has a notable effect on the vertical movement of water. In summer, there is also a significant difference at shallow depths within the permafrost region, indicating that simulations of deep-layer freeze-thaw processes play an important role in TSM in the shallow layers. In winter and autumn, the sensitivity of the runoff and groundwater scheme (RUN) gradually becomes dominant as depth increases because it also has an important effect on soil water movement (Figures 15c, 15d, 18c, and 18d). Furthermore, the energy-related STC scheme shows strong sensitivity in the northwest of permafrost regions in spring (Figures 16b–16d), which is again related to freeze-thaw processes; SFC determines the difference in evaporation through simulated energy differences, which have a strong sensitivity in summer and autumn in the arid Qaidam Basin. In addition, the SWC in spring and autumn is most sensitive to vegetation schemes in the regions covered with forests and grasslands (Figures 11b, 11c, 13a, and 13b), while the TSM is sensitive to INF and RUN schemes (Figures 16b, 16c, 18a, and 18b), which indicates these vegetation schemes have less effect on the soil ice simulations.

In regionally average simulations, INF, STC, RUN, and SFC schemes in winter; INF, STC, SFC, RUN, CRS, and DVEG in spring; INF, RUN, SFC, CRS, DVEG, and RAD in summer; and INF, RUN, SFC, DVEG, and CRS in autumn also have high sensitivity. Figure 19 (for the specific values of this sensitivity, see Table S3 in Supporting Information S1) shows the INF scheme had high sensitivity in TSM simulations throughout the year, especially in spring when the frozen soil begins to thaw. Meanwhile, the sensitivity of INF has a time lag between layers, such that the maximum value of sensitivity of the top soil occurs in the spring and is delayed to the summer as soil depth increases. Except for the INF and RUN schemes, the other schemes have lower sensitivity than other soil hydrothermal variables, which highlights the importance of these two schemes in soil water transport. In addition, the INF scheme has a great influence on the phase transition of soil water and ice because it has different sensitivities in SWC and TSM simulations.

Figure 19. Regionally average sensitivity of total soil moisture simulations in (a) winter, (b) spring, (c) summer, and (d) autumn. S is the Taylor Score.

From these results, it can be seen that the GT and ST simulations have high sensitivity to energy-related schemes across the TP. In simulations with dense vegetation, vegetation-related schemes also show high sensitivity. The simulations of SWC across the TP are not only controlled by schemes related to water movement but are also influenced by schemes related to energy transfer and vegetation variation because of the highly coupled freeze-thaw and evapotranspiration processes. Our results are similar to the sensitivity analyses of X. Li et al. (2021) and S. Sun et al. (2022), which were based on site observations. However, they concluded that INF has less effect in summer, which may be related to a different soil texture input within the simulations. When considering the total amount of solid and liquid water in soil, the effect of the freeze-thaw processes on the simulations is mainly reflected in the effect of vertical water transport. Meanwhile, different soil hydrothermal variables have different sensitivities.

Spatial Patterns of Suitable Parameterization Schemes

On the basis of remotely sensed MODIS land surface temperature and SMAP soil moisture data, we evaluated the applicability of schemes in the Noah-MP for simulating GT and shallow-layer SWC across the TP in different seasons. Noah-MP has many parameterization schemes that can be selected, and some schemes have temporal and spatial variations; however, some schemes are not sufficiently sensitive to create large differences in simulations, and it is difficult to describe each of them in detail. Therefore, only the most sensitive schemes outlined in Section 3.2 are discussed below; figures showing all the Noah-MP schemes and their proportions can be found in Figures S1–S8 in Supporting Information S1.

Suitable Options for GT Simulations

As described in Section 3.2.1 (Figures 5–8), the most sensitive schemes in the GT simulations of the TP are the energy-related SFC and ALB schemes and the vegetation-related DVEG and CRS schemes. Figure 20 displays a selection of different schemes based on the MODIS land surface temperature.

Figure 20. Selection of options from different schemes for ground temperature simulations using Moderate Resolution Imaging Spectroradiometer land surface temperature data based on the highest Taylor Score. (a–d) Are SFC; (e–h) are ALB; (i–l) are DVEG; (m–p) are CRS in winter, spring, summer, and autumn, respectively. Different colors correspond to different options.

The SFC is most sensitive across the western and northern TP throughout the year. However, different seasons require different options to match the remote sensing data. In winter, the Chen97 option is suitable for most areas except the central part of the TP (Figure 20a), while in spring, the M-O option is more suitable for the southern part of the TP (Figure 20b). In summer, when the vegetation growth is most vigorous, the M-O option is used across the entire TP except for the dry Qaidam Basin (Figure 20c). After the period of maximum vegetation growth, autumn sees the distribution characteristics return to a similar pattern as that seen in spring (Figure 20d).

The selection of the ALB scheme plays an important role in winter, and the CLASS option, which produces a smaller snow surface albedo due to its stronger aging effects, gives simulation results closest to the remote sensing data across the entire TP (Figure 20e).

Vegetation-related schemes (DVEG and CRS) have high sensitivity across the southeastern TP, and they extend to the central TP in summer and autumn. By comparing with remote sensing data, the corresponding options for DVEG and CRS in winter are 2 (on) and 1 (the Ball-Berry), respectively, except for the dry Qaidam Basin (Figures 20i and 20m). While the opposite options (DVEG is off and CRS is the Jarvis) should be applied in the southeastern TP with the variation of vegetation from spring to autumn (Figures 20j–20l and 20n–20p), the DVEG and the Ball-Berry options are more suitable for GT simulations in other areas of the TP. Notably, DVEG must be set to off when CRS is the Jarvis scheme based on the Noah-MP model's structure so that this opposite pattern of CRS and DVEG may appear when the Jarvis option is more applicable.

Suitable Options for Near Surface SWC Simulations

SMAP remote sensing can obtain the SWC of the shallow soil layer. We compared our simulation results of the upper soil layer (5 cm) with these SMAP data. Owing to insufficient remote sensing data in winter, only small areas in the southeast and north of the TP had enough data for model validation (Figure S5 in Supporting Information S1). Therefore, only the simulation results of the other three seasons are discussed herein.

Our simulation results show that the most sensitive schemes in the simulations of shallow-layer SWC across the TP include those related to energy transfer (STC and SFC), vegetation growth (DVEG), and freeze-thaw processes (INF and FRZ) (Figures 10–13).

Among these schemes, STC is most sensitive in most areas of the TP in winter and spring; exceptions are the southeastern and northwestern parts in winter and the southeastern part of the TP and the northern basin in spring. In spring (Figure 21b), seasonal freezing soil areas on the southeastern and northwestern TP are more suitable for option 1 (the semi-implicit scheme), while the southwestern and northern TP are more suitable for option 2 (the fully implicit scheme). SFC has a high sensitivity to the western TP in autumn, and the Chen97 option is recommended for these areas, but the recommended options in other regions are not consistent (Figure 21h).

Figure 21. Selection of options from different schemes for soil water content (at 5 cm) simulations using Soil Moisture Active Passive soil moisture data based on the highest Taylor Score. (a–d) are STC; (e–h) are SFC; (i–l) are DVEG; (m–p) are INF; (q–t) are FRZ in winter, spring, summer, and autumn, respectively. Different colors correspond to different options.

When vegetation is most sensitive in autumn, there is no uniform DVEG option across the southeast of the TP; instead, there is a mixed distribution of the two options (Figure 21l). Furthermore, in most areas of the TP during the summer and autumn when the vegetation is lush, option 2 (the DVEG model) can provide better simulation results than preset vegetation conditions.

In the simulations of surface layer SWC in summer, the scheme of frozen soil permeability (INF) is highly sensitive, but the two options are the same in the soil layers without ice, which indicates that the difference may come from the effect that the permeability of the deeper layers exerts on the shallower layers. Figures 21o and 21p shows that the Koren99 option is suitable in summer and autumn for permafrost regions. However, the NY06 option results in better performance in the northwestern permafrost regions in spring (Figure 21n). Many of the LSMs systematically underestimate top-layer SWC (Bi et al., 2016; K. Yang et al., 2009). The scheme with a greater impact on frozen soil permeability results in better simulations because it would reduce downward water transfer, making up for some of the underestimation. The supercooled liquid water scheme (FRZ) is most sensitive on the northern plateau in autumn, and the Koren99 option, which considers the increased interface between soil particles and liquid water, is most suitable for this season and in most regions of the TP (Figure 21t).

Our results show similar selections of options by comparing with the site observation results of X. Li et al. (2021), while some discrepancies still exist, which may be caused by uncertainty in forcing (Liang et al., 2019), soil texture and verification data (Dong & Ochsner, 2018; M. X. Li & Ma, 2010), and owing to the simulated and remote sensed variables having different characteristics under various land cover types (Brunsell et al., 2020). As seen from Figures S1–S8 in Supporting Information S1, the distribution characteristics of each of the better parameter options are not uniform, and there are differences between seasons and regions among different parameterized schemes. This means that it is difficult to accurately describe land surface processes that have spatiotemporal variations using only one scheme, as each scheme has its own advantages and limitations. Meanwhile, by comparing the results for GT and shallow layer SWC, we find that the optimal schemes for simulating different variables are not consistent; that is, schemes for improving GT simulation cannot guarantee the improvement of SWC simulation, and vice versa.

Conclusions

By utilizing Noah-MP multi-parameterization ensemble simulations, we examined the spatiotemporal patterns of sensitivities for simulating GT, ST, SWC, and TSM (including soil water and ice) in different soil layers across the TP. Furthermore, we revealed the relationships between parameterization schemes and climatic regimes, permafrost, and vegetation distribution across this region. Before the sensitivity analysis, we tested the spin-up time for soil variables across the TP, which allowed the model to reach an equilibrium state before simulations were performed. On this basis, suitable options for parameterization schemes for the entire TP region were selected by comparing the simulation results with remote sensing data. The main conclusions of this study are as follows:

• The spin-up times for GT, ST, SWC, and TSM (including soil water and ice) simulations are different. They show both horizontal (different regions) and vertical (different soil depths) differences. On the whole, the time required by ST is the longest, followed by SWC, and the time required by TSM is the shortest. The reason behind these time differences is that it takes a considerable amount of time for the soil freeze-thaw process to reach an equilibrium state; hence, a shorter time is required for TSM simulations that do not consider the phase change of water. Meanwhile, the spin-up times are longer in permafrost regions and shorter in both frozen and unfrozen regions.

• Sensitivity analysis of 12 different types of parameterization schemes in the Noah-MP LSM shows that the GT, ST, SWC, and TSM have different sensitivities, which in turn have notable temporal and spatial (including over different regions and different depths) variations. GT and ST are mainly sensitive to energy-related schemes, such as the snow and ST time scheme, surface layer drag coefficient, and lower boundary condition of ST. Vegetation-related schemes, such as DVEG and canopy stomatal resistance, play an important role after the vegetation growing season begins on the southeastern TP. In the sensitivity analysis of SWC, schemes related to both water (e.g., frozen soil permeability) and energy transport (e.g., snow and ST time scheme) show strong sensitivities, indicating that these schemes have substantial impacts on soil freeze-thaw processes and evaporation. However, the sensitivity of these energy-related schemes is weakened when simulating TSM (including the total amount of water and ice), indicating that the effect of direct phase change is reduced in these simulations so that they are only affected by the water transport schemes (e.g., frozen soil permeability) as a whole.

• Comparisons with MODIS land surface temperature and SMAP soil moisture data show that the selection of each parameterization scheme is not consistent over the entire TP across different seasons. With differences in climate regime, frozen soil, and vegetation distribution and variation, the selection of parameterization options also presents clear spatiotemporal differences. This means that different schemes each have their own advantages and disadvantages, but it is difficult for a single set of schemes to bring better performance with high applicability. At the same time, parameterization schemes that can improve GT simulation may not improve the simulation of shallow SWC, and vice versa.

This study focused on identifying and quantifying the sensitivities of parameterizations in the original Noah-MP. Owing to computational considerations, interactions between parameterization schemes were not explicitly addressed. Resolution of this issue will further advance our understanding of land surface processes and improve the performance of the model.

Acknowledgments

This research was funded by the Second Tibetan Plateau Scientific Expedition and Research (STEP) program (Grant 2019QZKK0103), National Program on Key Basic Research Project (2018YFC1505701), the Strategic Priority Research Program of Chinese Academy of Sciences (XDA20060101), and the National Natural Science Foundation of China (41830650, 91737205, 91837208, and 41905012). We thank the China Scholarship Council (CSC) for funding the author to study abroad, the Texas Advanced Computing Center (TACC) for providing us with computational resources and the reviewers for their valuable comments and experiences in refining the paper.

Data Availability Statement

The Noah-MP land surface model (version 4.1) is available at the website of the Research and Applications Laboratory at the National Center for Atmospheric Research (https://ral.ucar.edu/solutions/products/noah-multiparameterization-land-surface-model-noah-mp-lsm) and included in Niu et al. (2011). The China meteorological forcing data set is available at the National Tibetan Plateau Data Center (https://data.tpdc.ac.cn/en/data/8028b944-daaa-4511-8769-965612652c49) and stored in this in-text data citation reference: K. Yang and He (2019). The soil property data set is available at the National Tibetan Plateau Data Center (https://data.tpdc.ac.cn/en/data/11573187-fd64-47b1-81a6-0c7c224112a0) and included in Shangguan et al. (2013). The Moderate Resolution Imaging Spectroradiometer data set is available at (https://modis.gsfc.nasa.gov/data/dataprod/mod11.php), and the Soil Moisture Active Passive data set is available at (https://nsidc.org/data/smap/smap-data.html). The model results are listed at the National Tibetan Plateau Data Center (http://store.tpdc.ac.cn/index.php/s/ncJrCZWyPm7bDqb).