Bias Correction for ERA5-Land Soil Moisture Product Using Variational Mode Decomposition in the Permafrost Region of Qinghai–Tibet Plateau

Soil moisture (SM) is one of the key measures to understand the land-atmosphere interaction and permafrost dynamics in the Qinghai–Tibet Plateau (QTP). ERA5-Land is a new reanalysis product with high spatial resolution (9 km), which can provide long-term SM data with a large spatial coverage as well as at multilayer soil depths. However, preliminary comparisons with in-situ data show that the ERA5-Land SM product generally underestimates the seasonal variability and demonstrate a positive bias on the QTP. In this article, we proposed to utilize the mode decomposition method to correct such bias. Specifically, through using the variational mode decomposition, we decomposed the long time-series of ERA5-Land SM data into a series of intrinsic mode functions, and found that the SM seasonal variation can be well represented by the low-frequency modes, which were then selected to feed a regression model for the bias correction. The single-site bias correction results show that our method significantly improves the accuracy of ERA5-Land SM product with bias reduced by 0.22 m3/m3, 0.3 m3/m3, and 0.15 m3/m3 for alpine meadow, alpine steppe, and alpine desert sites, respectively. Together with the slightly reduced accuracy but still acceptable results for the cross-site bias correction, we successfully demonstrate the potential of the mode decomposition method for the bias correction of the ERA5-Land SM product at regional scale. Our method is of great use to study climate impact on regional ecohydrological processes and the permafrost changes in the QTP region.


I. INTRODUCTION
T HE Qinghai-Tibetan Plateau (QTP), often known as the Third Pole, contains the world's largest permafrost area in the middle and low latitudes, which accounts for ∼40% of the total area of QTP [1]. Widespread permafrost degradation on the QTP has been observed in pace with the increased air temperature from 1980 to 2018, where the active layer has thickened at a rate of 19.5 cm per decade [2]. The annual freeze and thaw of the active layer controls the soil water movement and redistribution [3], [4]. Therefore, increasing thaw depth and permafrost degradation will have significant impact on surfacesubsurface water interaction and runoff production, thus induce regional soil moisture (SM) variations, and alter the groundwater recharge [5].
On the other hand, SM, as a vital component of the terrestrial water cycle, regulates the exchange of water and energy between the atmosphere and the land surface [6]. SM also affects the soil hydraulic and thermal properties, and serves as one of the critical factors impacting the soil freezing-thawing process [7]. SM changes have a direct impact on the ecosystem processes in the permafrost regions, affecting the decomposition and accumulation of soil carbon, greenhouse gas emission, and the living environment of vegetation [8], [9]. Therefore, long-term and reliable SM data are useful for the understanding of the land-atmosphere interaction and the environmental evolution in the permafrost region.
Due to the harsh environment, long-term SM field observation networks are extremely difficult to maintain on the QTP. Remote sensing retrievals [10], land hydrological modeling [11], and reanalysis data [12] can provide continuous information of large-scale SM variations. However, remote sensing technology only provides near-surface SM information with the sensing depth generally less than 10 cm [13]. On the other hand, the accuracy of land hydrological models is limited by deficiencies in the model structure, model parameterization scheme, and uncertainties in the meteorological forcing data [14]. In contrast, the reanalysis products provide comprehensive data that integrate ground observations, multisource satellite remote sensing data, and numerical model simulations [15]. The advantages of the reanalysis products have been demonstrated for its capability in capturing spatial and temporal SM changes across global [16]. This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ ERA5-Land is a state-of-the-art reanalysis product from the European Centre for Medium-Range Weather Forecasts (ECMWF), with a high spatial resolution (9 km), which provides SM product at four layers (0-7, 7-28, 28-100, and 100-289 cm) [17]. Some studies have preliminarily assessed the performance of the ERA5-Land SM product. The correlation is relatively high between the ERA5-Land surface (0-7 cm) SM and 826 in-situ measurements, primarily located in the USA and Europe, with a median correlation (R) of 0.72 [18]. On the QTP, the ERA5-Land product captures the seasonal fluctuations in SM, but underestimates the magnitude of temporal changes in SM [19]. Additional comparisons indicate that the ERA5-Land product generally overestimates SM at 0-100 cm in relation to 16 in-situ sites in the permafrost region, with a median bias of 0.17 m 3 /m 3 and a median unbiased root mean square error (ubRMSE) of 0.03 m 3 /m 3 from July to August during 2016-2018 [20]. The low ubRMSE value implies that the ERA5-Land product performs better if the bias is removed. In order to obtain reliable SM data to support the study of permafrost changes, it is necessary to remove the bias in the current SM products prior to any further applications.
Commonly used bias correction methods include linear regression analysis [21], mean and standard deviation (μ-σ) matching [22], min-max correction [23], cumulative distribution function (CDF) matching [24], etc. Linear regression is a simple and effective method, but it requires a good linear correlation between the original data and the reference data [25]. Min-max correction and μ-σ matching cannot improve the correlations between the two datasets [26]. The CDF estimation requires a long data record and may ignore transient or seasonal dynamics [27]. In order to overcome the above-mentioned issues, signal decomposition methods have been also introduced for bias correction. This type of technique can decompose the signal into multifrequency (or resolution) forms that are easy to analyze and characterize, and consequently separate redundant components or interference signals from the original signal [28]. Haerter et al. [29] proposed to decompose the time series of atmospheric variables in order to independently correct the errors of climate model outputs at different time scales. Kusumastuti et al. [30] used discrete wavelet transform for bias correction of variability and trends in climate model simulations. In the realm of signal decomposition, the empirical mode decomposition (EMD) and its variations form a class of adaptive methods that are suitable for nonlinear and nonstationary hydrological time series [31]. Variational mode decomposition (VMD) is a novel mode decomposition method, which is based on a strict mathematical foundation, overcoming the problems of end-point effect, mode-mixing phenomenon, and noise sensitivity in the class of EMD methods [32], [33]. The application of VMD in the hydrological data series mainly includes river runoff prediction and trend analysis [34], [35].
In this study, we propose applying the VMD method to the bias correction of the ERA5-Land SM product in order to obtain more reliable SM data at multilayer depths and during a multiyear period in the QTP permafrost region. Our goal is to extract the key components of the SM signals using VMD, to restore the seasonal signals of SM variability, and to eventually remove the bias through a regression model. The rest of this article is organized as follows. In Section II, we describe the data used in this study. In Section III, our processing procedure is introduced in detail following a basic description of VMD. In Section IV, we report the results of single-site and cross-site bias corrections that were performed for three main vegetation types on the QTP, including alpine meadow, alpine steppe, and alpine desert, using the correction equations obtained from the site itself and other sites with the same vegetation type, respectively. We discuss the results in Section V, followed by the conclusion.

A. In-Situ Measurements
The permafrost area of QTP is mainly distributed in the areas with the elevation range of 4000∼5500 m a.s.l, with the mean annual temperature lower than −2.5°and the mean active layer thickness of ∼2.3 m [36], [37]. The precipitation is mainly concentrated during the monsoon season, with an average annual precipitation ranging from 99 to 600 mm [38]. The dominant land cover types ( Fig. S1 in the supplementary material) include alpine swamp meadow, alpine meadow, alpine steppe, alpine desert, and bare land [39].
Our study area is mostly located in the QTP hinterland permafrost region. The in-situ observation data used in this article mainly cover three vegetation types (alpine meadow, alpine steppe, and alpine desert), as illustrated in Fig. 1. The data at D66, D105, and D110 sites are from the SM and soil temperature network of the coordinated enhanced observing period Asia-Australia monsoon project (CAMP) in the Tibetan Plateau (CAMP/Tibet, 2001-2005), whose main task was to study the water cycle over the Tibetan Plateau [40]. The data from other sites are provided by the Chinese Academy of Sciences the Cryosphere Research Station (CRS) on the QTP, Chinese Academy of Sciences, which has been developing a soil temperature and moisture monitoring network over the permafrost region of QTP since 2004 [41]. SM was measured using the Steven Hydra Probe (Cambell, USA) at multiple depths and a 30-min time interval with an accuracy of ±2.5% [41]. More details about the in-situ sites are given in Table I.

B. ERA5-Land Reanalysis Product
The latest ERA5-Land product is an enhanced global reanalysis dataset for the land surface, being produced by the ECMWF within the framework of the Copernicus Climate Change Service (C3S) of the European Commission [17]. The product is obtained by the global high-resolution numerical integrations of the Hydrology Tiled ECMWF Scheme for Surface Exchanges over Land (HTESSEL), which is driven by meteorological forcing downscaled from the ERA5 reanalysis, including elevation correction for thermal near-surface states [17]. In the soil hydrology scheme [42], the soil heat budget follows the Fourier diffusion law. The soil volumetric heat capacity is set as a constant and the heat conductivity is determined by the total soil water content, with the conductivity changes caused by soil water phase change ignored. Precipitation is collected in the interceptor reservoir until saturation, and excess precipitation is divided into surface runoff and infiltration. The soil discretization uses a four-layer soil system, overlain by a single layer of snow. The vertical movement of soil water is described using Richards' equation [43]. The top boundary condition is determined by infiltration and surface evaporation, and free drainage is assumed at the bottom. ERA5-Land provides hourly data from 1979 to the present with the highest spatial resolution (∼9 km) among current reanalysis products. The data are available from 1 . In this study, we use the volumetric water content (VWC) data for the four layers of soil column (0-7, 7-28, 28-100, and 100-289 cm).

C. Ancillary Datasets
Additional datasets used in this article include 1-km vegetation map of QTP (Fig. S1) [39] and the EVI (enhanced vegetation index, MCD43A4) during the growing season (June to September). The vegetation types of the measurement sites were extracted from the vegetation map, which cover three main vegetation types in the QTP permafrost region [39]. In this study, the mean EVI during the growing season from 2000 to 2020 was used to represent the vegetation coverage. The EVI map with a spatial resolution of 250 m on the QTP was generated from the MODIS 16-day composite EVI dataset (MOD13Q1) [44] using Google Earth Engine. We reclassified the vegetation type of a site if it shows a large discrepancy in terms of the growing-season EVI, from other sites with similar vegetation types. For example, QT08 is originally classified as an alpine meadow site in the vegetation map. However, the mean growing-season EVI of this site is 0.10, which is much smaller than other alpine meadow sites (mean growing-season EVI of 0.22). Since the growing-season EVI at alpine steppe sites ranges from 0.09 to 0.16, the QT08 site is reclassified as alpine steppe.

III. METHODS
In this study, we propose to utilize the mode decomposition method to correct the bias in the ERA5-Land SM data. In this section, we take the VMD, one of the most popular mode decomposition methods, as an example to describe the methodology.

A. Variational Mode Decomposition Method
VMD is a new, intrinsic and adaptive variational signal decomposition technique proposed by Dragomiretskiy and Zosso [32] to decompose an input signal f (t) into K discrete number sets of band-limited intrinsic mode functions (IMFs) with specific sparsity properties nonrecursively. The constrained variational problem described by VMD for the f (t) time series is where K is the number of modes to be decomposed, mode compacts around a corresponding center pulsation ω k , i.e., {ω k } = {ω 1 , ω 2 , . . . , ω k } . This center frequency can be updated in the Fourier domain, representing the frequency response of the decomposed subsignals. δ(t) is the Dirac distribution, and * represents the convolution operation. The quadratic penalty term α and Lagrange multiplier λ are introduced to find the optimal solution, which transforms the constrained problem into an unconstrained one. The quadratic penalty allows the reconstruction with high fidelity, typically in the presence of additive Gaussian noise. Therefore, the augmented Lagrange is introduced as follows: The decomposed modes and center frequencies are solved by the alternate direction method of multipliers [45]. A detailed description and solution of VMD are provided by Dragomiretskiy and Zosso [32].

B. Bias Correction Procedure and Evaluation Metrics
The specific procedure for the SM bias correction using VMD is depicted in Fig. 2. Our goal is to match the seasonal amplitude of ERA5-Land SM data with the in-situ measurements while reducing the overall bias. The in-situ SM data in the permafrost region are generally associated with large errors due to logistics challenges [46]. In addition, the sensor may be affected by the soil freeze-thaw cycle, which would result in sensor position drift and malfunction, affecting the accuracy of the measurements [47]. Therefore, we first performed data screening of the in-situ observations. Specifically, we excluded the sites with obvious errors, such as the sites or data with almost identical SM values at two different depths. We then averaged the hourly data of CAMP/Tibet and ERA5-Land data to the daily time step as CRS SM data. We adopted the pixel-to-point matching strategy used by many studies in data-scarce regions [48], through selecting the ERA5-Land grid nearest to the in-situ site with a distance ranging from 2.55 to 6.35 km. To match the corresponding soil layer between ERA5-Land and multilayer data, we first applied the linear interpolation to the multidepth in-situ SM data to obtain the SM profile between the uppermost and deepest measurement depths, with vertical resolution of 1 cm. Based on the SM profile, we calculated the mean SM of the soil layer corresponding to the ERA5-Land layer depth as follows: where SM l 1 −l n is the average SM in the soil layer of l 1 to l n cm (m 3 /m 3 ); l 1 , l n refer to the upper and lower soil depth (cm) of a specific ERA5-Land soil layer, respectively. SM i refers to the interpolated in-situ SM at the i th depth (m 3 /m 3 ); u, d are the uppermost and deepest depth of the in-situ SM profile within the depth range of l 1 − l n (cm). If the SM profile completely covers the specific soil layer of ERA5-Land, u = l 1 , d = l n . Then, the difference time series, D(t) , was obtained by subtracting the ERA5-Land SM data from the in-situ SM data in the same depth range, as follows: where SM insitu (t) and SM ERA (t) are the in-situ SM and ERA5-Land SM, respectively. VMD is then used to decompose the difference sequence D(t) into K discrete numbers of subsignal modes u k (t), as follows: The disparities between the in-situ and ERA5-Land SM time series are difficult to assess directly due to the high variability in SM data at the daily scale. Such deviations can be represented by combining the suitable modes decomposed by VMD. In the current study, we set K equals 5. Our experiments show that regardless of this number of decomposed modes, the last two low-frequency modes generally have the largest amplitudes, indicating that the dominant differences between the ERA5-Land and in-situ data generally occur at the seasonal and annual time scale. On the contrary, the magnitudes of the decomposed high-frequency subsignals in the first three modes are relatively small and it is difficult to develop an empirical function to fit these highly nonlinear data, thus neglected in this study.
Therefore, u 4 (t) and u 5 (t) were selected for the regression analysis to obtain the correction equation. The modes u 4 (t) and u 5 (t) represent the bias at the seasonal and annual time scale, respectively. In this study, the sinusoidal and linear functions are implemented to fit u 4 (t) and u 5 (t) as in (6) and (7), respectively. Finally, the combined equation is applied to the ERA5-Land SM for bias correction, as in (8).
In this study, we performed both single-site and cross-site bias correction. For the single-site correction, we chose sites with a time span of more than three years. A regression equation was then developed based on ∼60% of the time series at this site. The remainder of the time series at the same site is then used for validation. As for the cross-site bias correction, the correction equation for the site is obtained from all other sites with the same vegetation type excluding that site. The purpose of the cross-site correction is to verify whether the bias characteristics of ERA5-Land SM are consistent among the sites with the same vegetation type, which is critical for regional SM bias correction.
In this study, the mean bias (Bias), correlation coefficient (R), and root mean square error (RMSE) are used to evaluate the accuracy of the bias correction at daily scale. The expressions of the three indicators are as follows: where N is the total number of samples; θ n andθ n are measured and estimated SM values, respectively.θ andθ are the corresponding mean values during the period.  Fig. 4(a). The overall fitting is good, with R ≥ 0.94, indicating that the equations selected for IMF4 and IMF5 are appropriate. Fig. 4(b) displays the absolute error distribution of initial ERA5-Land SM and corrected ERA5-Land SM at QT03. The error of the original ERA5-Land SM data averaged over the four layers is 0.19 m 3 /m 3 . The corrected SM data show much reduced error, with an average value of 0.04 m 3 /m 3 .

A. Bias Correction for Single Site
Results for all single-site bias correction for alpine meadow, steppe, and desert sites are shown in Figs. 5, 6, and 7, respectively. The evaluation statistics of the initial and corrected ERA5-Land SM are given in Table II. There is an overall good consistency between the corrected ERA5-Land SM and in-situ SM for all three vegetation types. Improvements in the R, Bias, and RMSE are seen for almost all the sites. For alpine meadow sites, the mean R increases from 0.66 to 0.93, while the bias and RMSE decrease from 0. The bias correction does not work well at the 7-28 cm and 100-289 cm layers of ZNH site, as well as the 100-289 cm layer of TSH site (alpine desert). The dynamics of the in-situ SM time series during summer at these two sites are quite different from other sites, which might be due to lateral water flow from nearby lakes [49]. At the ZNH site, the in-situ data show that SM increases significantly from August to October (see Fig. 6), an increase of ∼250% to prior SM at each soil layer. ZNH site is located near to the Zonag Lake, and the SM may be significantly impacted by the lake runoff; the SM may be also influenced by the outburst event of Zonag Lake which altered the regional hydraulic connectivity and accelerated permafrost degradation [50]. As shown in Fig. 7, the SM at layer of 100-289 cm at the TSH site increases dramatically from July to December every year, nearly four times of the SM at layer of 28-100 cm. Therefore, the high SM at layer of 100-289 cm should not be entirely supplied from infiltration of soil water in the upper layer and may be caused by the lateral flow from the neighboring Tianshuihai Lake, about 80 m away from the site. Through comparing the soil temperature and SM profiles of the TSH site (Fig. S2(a) and (b) in the supplementary material), we found that SM is unusually high even when soil temperature drops below 0°C at the soil depths below 150 cm. The unfrozen water curve (Fig. S2(c) in the supplementary material) indicates that the freezing point at depths of 160 and 180 cm at the site is below −1°C. Therefore, when soil temperature is above or very close to the freezing point, the unfrozen water remains high. The presence of salt, capillary action from fine soil pores, or high pressure can cause the freezing point depression [51]. Strong upward shortwave radiation is present at the TSH site year-round due to the dry and salt-rich ground surface caused by large evaporation and low precipitation [41]. During the thawing process of the active layer, SM and salt mainly migrate downward, so the salt content in the deep soil layer increases [52]. Therefore, high salt content at depths of 160 and 180 cm may result in the freezing point depression at the TSH site.  Table III for all the sites with different vegetation types. The results show that after the correction, the bias is significantly reduced for all types of land covers. The mean R, bias, and RMSE of alpine meadow sites are 0.92, 0.03 m 3 /m 3 , and 0.08 m 3 /m 3 , respectively. For alpine steppe sites, the R of corrected SM ranges from 0.58 to 0.95, with a mean value of 0.84. The average bias and RMSE are 0.07 m 3 /m 3 and 0.09 m 3 /m 3 , respectively. The biases for all sites except QT05 are relatively small, i.e., less than 0.06 m 3 /m 3 . The R is increased from 0.04 to 0.86 after the correction at alpine desert sites. The average values of bias and RMSE at alpine desert sites are 0.00 m 3 /m 3 and 0.03 m 3 /m 3 , respectively. Note that the accuracy of measured SM is ±0.025 m 3 /m 3 [41], which is close to the mean RMSE of our corrected data. Moreover, the R of the SM times series increases for all sites after bias correction, so the cross-site bias correction results are considered acceptable.

B. Cross-Site Bias Correction and Validation
However, cross-site bias correction may not work well at sites show different SM variability from other sites that were used to develop the correction equation. For example, the in-situ SM at QT01 shows a seasonal decaying trend, quite different from other alpine meadow sites (Fig. 5), and the errors of original ERA5-Land SM at QT05 are much larger than other alpine steppe sites (Table II). Although the corrected SM and in-situ SM have a strong correlation and the errors are reduced after correction for the QT01 and QT05 sites, the bias and RMSE still remain high (Table III). It indicates that the inconsistency of SM variations and the discrepancy of ERA5-Land biases from other sites with the same vegetation type may be the main reason for the large errors of cross-site correction at these two sites. Therefore, these two sites are not included in the cross-validation of alpine meadow and alpine steppe types. We compared the error statistics for the sites included in crosssite and single-site corrections, as shown in Fig. 9. It shows that the accuracy of single-site correction is better than the cross-site bias correction for both alpine meadow and alpine steppe sites. The average R between cross-site corrected SM and in-situ SM is slightly lower than that of single-site correction. The mean bias and RMSE of the cross-site correction for the alpine meadow are higher than that of the single-site correction by 0.03 m 3 /m 3 and 0.04 m 3 /m 3 , respectively. For alpine steppe sites, compared with single-site correction results, the average bias and RMSE of cross-site correction are increased by 0.05 m 3 /m 3 and 0.06 m 3 /m 3 , respectively. Moreover, the single-site validation shows that the proposed method performs well in both alpine meadow and alpine steppe sites. However, in the cross-site validation, our method works better at the alpine meadow sites than at the alpine steppe sites. This is due to the higher correlation coefficient and smaller biases between the original ERA5-Land and in-situ SM data at alpine meadow sites than the alpine sites, which indicates that ERA5-Land performs relatively well in humid regions, but may have larger uncertainties in arid regions.

V. DISCUSSION
The bias corrected ERA5-Land SM data generally provide more accurate information on the temporal pattern of SM than the initial ERA5-Land SM data. The bias corrected SM data exhibit a stronger correlation and lower errors with the in-situ data than the initial ERA5-Land SM data at the majority of sites. The most significant feature of our method is that it corrects the seasonal change of SM while reducing the overall bias. In addition, our method is applicable to both surface and deep SM. The cross-site validation results show similar accuracy as the single-site bias correction method, with slightly better performance at alpine meadow sites than at alpine steppe sites, which indicates that our method's applicability at regional scale. However, there are also some limitations in our method as discussed in the following.
Ignoring the high-frequency signals (IMF1-3) in the SM differences between the ERA5-Land and in-situ SM data reduces the accuracy of the corrected SM data to some extent. We only selected the low-frequency modes, i.e., IMF4 and IMF5, for bias correction. However, the high-frequency modes may also have a considerable contribution to the total bias, such as at the Ch06 site. The high-frequency signal mainly reflects the dynamic characteristics of SM with high temporal resolution (e.g., at daily scale). SM is highly dynamic over time, especially at the shallow surface, which is susceptible to precipitation inputs and changes in surface evaporation. Therefore, ignoring the high-frequency signals may also lead to errors more pronounced in the surface layer than that in the deep layer, and shall be addressed in the future studies.
Using a simple sinusoidal function (6) to recover the seasonal amplitude of ERA5-Land SM may result in overestimation of unfrozen water content during soil freezing. The corrected SM is controlled by the sinusoidal function and generally shows a "V" shape, while the SM in permafrost regions is affected by the soil freeze-thaw process and the in-situ SM series generally show a "U" shape. This is particularly true for the surface soils due to its rapid freezing process, and sudden drops in the unfrozen water content when the soil temperature drops below the freezing point [53]. The SM at deeper soils, however, may show a gradual decline due to the latent heat release during the freezing process, which can drag soil temperature close to the freezing point for a long period of time (i.e., the zero-curtain period) [54]. As a whole, the unfrozen water changes during these processes may not be well represented by the sinusoidal function. In addition, the ERA5-Land does not distinguish liquid water and ice in the SM product, while the in-situ data only represent the soil liquid water. This may explain why the sinusoidal function may not represent well the freeze-thaw effects on the SM changes. To evaluate the performance of our method during the freezing period, we calculated the errors during this period. The first day when the soil temperature drops below or above 0°for 3 consecutive days is defined as the onset date of soil freezing or thawing, respectively [55]. The average RMSEs for alpine meadow, alpine steppe, and alpine desert sites during the freezing period are 0.03 m 3 /m 3 , 0.02 m 3 /m 3 , and 0.03 m 3 /m 3 , accounting for 78%, 60%, and 72% of the errors during the entire season, respectively (see Table SI in the supplementary material). Therefore, the inability of the model to capture the characteristics of unfrozen water changes during the freeze-thaw period is one of the main error sources in our correction results. It is noted that the migration and changes of soil liquid water are important to understand the process of soil heat transfer, frost heave, and thermal-hydro-mechanical interaction [56]. Future work should include estimating the unfrozen water content from the ERA5-Land SM product and explore different types of functions for the SM bias correction.
Discontinuity in the variations of SM time series also introduces additional errors to the bias corrected SM data. In the single-site correction, because the time-series of SM at the site is relatively stable, our method works rather well through correcting the interannual bias decomposed from VMD. However, due to the inconsistencies of SM variations among different sites, the interannual bias of SM may not be well accounted using the cross-site bias correction method. As shown in Fig. S3, the in-situ SM data at QT01 for the shallow layer of 0-28 cm shows that the soil begins to dry up after 2005, with SM dropping by roughly 45% during summer. However, there are no significant decreases in SM at other meadow sites and ERA5-Land also fails to capture the interannual trend of SM at the QT01 site. As a result, the cross bias corrected ERA5-Land product overestimates SM at the QT01 site after 2005.
The performance of our bias correction method may also be hampered by the deficiencies of ERA5-Land. The accuracy of the reanalysis products depends on a mix of observations and model forecasts. In data-sparse regions, such as the QTP, the quality of the reanalysis data is strongly dependent on the quality of forcing data, model structure, and assumptions [57]. Some studies have shown that meteorological forcing had a more significant impact on surface SM simulations, while the LSMs play a more important role in the SM simulations at deeper depths [58]. Due to the inability of HTESSEL to represent horizontal soil water transport within soil layers [42], ERA5-Land is unable to capture the changes in the SM caused by lake inflow or other lateral water recharge at ZNH and TSH sites. Such SM changes may also result from the outburst of Zonag Lake in 2011, which changed the regional hydraulic connectivity and induced permafrost degradation [50]. On the other hand, LSMs employed in the global reanalysis products in general are not specifically designed for permafrost regions. It is still challenging for most LSMs to accurately describe the glacier and snow-melt, soil freezing and thawing, and other physical processes in the Cryosphere [36], [59], [60]. Therefore, considerable uncertainties may exist in the reanalysis product for those remote regions. For example, the ratio of soil ice content to unfrozen water content will constantly change at different stages of the freezing-thawing process, which significantly changes the soil thermodynamic and hydraulic properties [61]. However, the HTESSEL model ignores the influence of soil water phase changes on the soil thermal properties, which may lead to large deviation of soil temperature and unfrozen water simulations during the freezing period. In addition, HTESSEL assumes free drainage at the bottom of soil column, so it cannot well represent the restriction effects of underlying permafrost on the subsurface flow. Therefore, SM may have larger uncertainties at deeper soils in the permafrost region. This might explain the large bias of ERA5-Land SM product at the QT01 site. Due to the smaller active layer thickness of QT01, SM at deeper layers of QT01 should be mainly frozen and much smaller than that of other alpine meadow sites; however, ERA5-Land SM at deep soil layers is not consistent with the in-situ data at the QT01 site.   Moreover, the input parameters of the model may not be able to accurately distinguish the permafrost and seasonal frozen ground, so that the land model may have large errors in the transitional zone, such as at the QT05 site.
Our method may also face difficulties in the alpine desert areas where SM variability is much stronger than in alpine meadow and steppe areas. Fig. 10 shows the difference signal D(t) for the surface layer after VMD at the D66 site. The IMF4 does not conform to our expectations and is unable to match the sinusoidal equation that we specified, which indicates the variability of surface SM at the alpine desert site is so large that the ERA5-Land data could not effectively capture the seasonal changes of SM. This may be because that the meteorological forcing data used by ERA5-Land, as well as the model parameters related to soil hydrothermal transport, are not accurate enough to represent the hydrothermal conditions and the soil characteristics in the QTP arid area [62], [63]. ERA5-Land tends to overestimate the precipitation and underestimate the soil temperature on the QTP [48], [63], which affects the infiltration and evaporation process, resulting in SM overestimation during the thawing season. Soil hydrological parameters in the HTESSEL model are determined by the Food and Agriculture Organization (FAO) dataset [42], which tends to overestimate the clay content on the QTP, likely resulting in higher SM estimates [6], [62].
In addition to the above-mentioned issues, due to the sparse in-situ data and limited depth information, the scale mismatch in both spatial and vertical representation between the in-situ data and ERA5-Land SM product may also introduce uncertainties to our results. The variation of SM with depth in the permafrost region is generally not linear [64], so the use of linear interpolation to obtain the SM profile likely adds additional uncertainty. Moreover, almost all the in-situ data from the alpine meadow and alpine steppe sites were collected on the relatively flat terrain along the Qinghai-Tibet Railway. Therefore, the applicability of our method in other regions has not yet been verified, but should be included in the future with additional available in-situ observations.

VI. CONCLUSION
In this study, we have proposed to utilize the mode decomposition methods, and more specifically, have implemented VMD to correct the bias in the ERA5-Land SM product relative to the in-situ SM data in the QTP permafrost region. Through comprehensive single-site validations, our bias correction method is demonstrated to be able to extract the deviations of ERA5-Land SM products for most sites. The temporal correlation between the ERA5-Land SM and in-situ SM is significantly improved with overall bias reduced that generally exists in the ERA5-Land SM data. The results of cross-site bias correction also demonstrate an improvement in the consistency between the ERA5-Land and in-situ SM data, indicating this method's potential for regional-scale SM correction. The ERA5-Land SM after the cross-site correction at the alpine meadow sites generally has a higher accuracy than that at the alpine steppe sites, which may indicate possible deficiency of ERA5-Land model in the more arid regions. Therefore, developing improved LSMs and SM retrieval algorithms suitable for QTP permafrost regions are essential for accurate SM estimates. Our method can be applied to various remote sensing and reanalysis products to extend the in-situ SM data with limited time length or large gaps, for analyzing the characteristics and variations of SM at local scale, and even for supporting large-scale studies of hydrothermal characteristics of QTP permafrost.

AUTHOR CONTRIBUTIONS
Rongxing Li led the project. Tian Chang and Yuliang Wen carried out the experiments. Tian Chang wrote the draft. Yonghong Yi supervised the study, interpreted the results, and wrote and edited this article. Tong Hao proposed and designed the methodology, analyzed the results, and finalized this article. All coauthors contributed to data interpretation and editing of this article.
Tian Chang received the B.S. degree in remote