On the Spin‐Up Strategy for Spatial Modeling of Permafrost Dynamics: A Case Study on the Qinghai‐Tibet Plateau

Spin‐up is essential to provide initial conditions for land surface models (LSM) when they cannot be given reliably as in the application to regional permafrost change studies. In this study, the impacts of spin‐up strategy including total spin‐up length and cycling scheme on modeling of permafrost dynamics on the Qinghai‐Tibet Plateau (QTP) were evaluated through two groups of experiments using a modified Noah LSM. The first group aims to test different total spin‐up lengths and the second group for different cycling schemes. The results show that the presence of permafrost prolongs the convergence of the model. Vertically, the slowest convergence is observed at the permafrost table. The insufficiency of total spin‐up length is prone to underestimate permafrost area and overestimate the degradation rate. Different cycling schemes considerably affect the resulting initial thermal fields and result in degradation rates with a difference of 3.37 × 103 km2/a on the QTP, which exceeds the difference (2.92 × 103 km2/a) in the degradation rates reported in existing studies. The multi‐year cycling scheme is generally preferred, but overlong cycle length should be avoided to prevent the introduction of climate change trends in the spin‐up period. We recommend a spin‐up strategy of a 500‐year cycling with the first 5‐ to 10‐year of forcing for modeling permafrost on the QTP with the Noah LSM. Our findings highlight the importance of the spin‐up strategy, which is usually neglected in present LSM‐based permafrost modeling studies.

. Numerous observation-based evidences show that permafrost has changed globally to varying degrees in recent decades (Biskaborn et al., 2019;Smith & Riseborough, 2002;Zhao et al., 2010). Modeling is a powerful tool for reconstructing, assessing and predicting permafrost changes in the past, present and future Walvoord & Kurylyk, 2016). In recent years, thanks to their explicit physics basis, land surface models (LSMs) have been continuously applied and improved for modeling thermal and hydrological processes in permafrost regions (Ekici et al., 2014;Lawrence & Slater, 2005; G. Zhang et al., 2019).
In an LSM, water and heat transport processes through the soil are typically described using diffusion equations (e.g., Lawrence et al., 2011;Niu et al., 2011). When we solve those nonlinear partial differential equations by a finite-difference method, the initial values of state variables, such as soil temperature and moisture content, must be specified in advance (Gubler et al., 2013). At site scale, accurate initial fields and model parameters can often be determined by in-situ measurements (Gao et al., 2015;Nicolsky et al., 2007). Lumped model parameters, or parameters that are hard to measure directly, can also be calibrated using long-term observations (Ekici et al., 2015;Westermann et al., 2016). In addition, many site-scale studies have been conducted to address the shortcomings of an LSM in simulating permafrost dynamics, such as improving the parameterization of ground freeze-thaw processes (Ekici et al., 2015;Sun et al., 2020), snow cover and snowpacks Li et al., 2021), and soil organic properties (Jafarov & Schaefer, 2016;S. Yang et al., 2021;Zhu et al., 2019). Therefore, the impact of initial conditions on site-specific simulation is much reduced (Gubler et al., 2013), although some hydrogeological studies in cold regions have also cautioned about the potential impact of model initial fields on groundwater contribution to streamflow (Bense et al., 2012) and the ground thawing process (Mckenzie & Voss, 2013). However, the initial fields often cannot be reliably provided in a spatial modeling application since observations at most locations are absent (Lawrence et al., 2008;Riseborough et al., 2008). Thus, spin-up becomes a plausible way to provide the initial fields for LSM simulations (Rodell et al., 2005; Z.-L. Yang et al., 1995).
After a sufficient spin-up procedure, the model should be internally self-consistent, that is, without artificial drift in model stores or in prognostic variables (Cosgrove et al., 2003;Ednie et al., 2008;Lamontagne-Hallé et al., 2020). At this point, the state variables reach quasi-equilibrium with the atmospheric conditions of the spin-up period. Many factors can affect a sufficient spin-up, such as the initial states provided to start the spin-up (Cosgrove et al., 2003;Shrestha & Houser, 2010), the numerical model used for the spin-up (Gao et al., 2015; Z.-L. Yang et al., 1995), and climatic characteristics of the study area (de Goncalves et al., 2006). Cold regions usually require a much longer total time length (e.g., hundreds of years) for spin-up than other regions de Goncalves et al., 2006), mainly due to repeated freeze-thaw cycles (Gao et al., 2015) and large soil depth involved (Lawrence et al., 2008). In permafrost regions, due to the "heat reservoir" effect of deep soils Sapriza-Azuri et al., 2018;Stevens et al., 2007), the modeling of permafrost require a large simulation depth, which is usually greater than the depth of zero annual amplitude (DZAA) of ground temperature, which averages over 10 m on the QTP (Xie et al., 2015). It increases the memories of soil temperature and moisture conditions and thus poses a challenge for a suitable spin-up strategy to obtain acceptable initial fields in permafrost regions Riseborough et al., 2008).
Given the disequilibrium between permafrost and current climate, a transient spin-up simulation is implemented by a regional climate model that includes an LSM  or by reconstructing long-term climate data dating back to hundreds of years ago (Ednie et al., 2008;Sapriza-Azuri et al., 2018). However, this method appears difficult when providing initial fields for regional permafrost simulations. The climate model operates on a coarse spatial resolution and is generally subject to great uncertainties (Slater & Lawrence, 2013;W. Wang, Rinke, et al., 2016). It is also unlikely that the reconstructed data will be available for regional coverage due to the lack of regional dendroclimatic information, for example, . As a compromise, the spin-up simulation can be driven by steady-state boundary conditions (Bense et al., 2012), that is, conditions do not vary with time and correspond to an average meteorological condition. A steady-state spin-up is limited in its assumption that the ground states at a study site are initially in equilibrium with the constant hydroclimatic conditions (Lamontagne-Hallé et al., 2020). If seasonality is considered, a dynamic equilibrium spin-up simulation (i.e., driven with seasonally varying meteorological forcing without any yearly changes) is preferred for permafrost modeling (Lawrence et al., 2008;Lee et al., 2014;G. Zhang, Nan, Yin, et al., 2021). However, for larger regions where conditions vary spatially, it becomes even more difficult to select a year with representative atmospheric conditions for both steady-state and dynamic equilibrium spin-up approaches.
In practice, we often constitute spin-up forcing of a certain length (i.e., the total spin-up length) by cycling the meteorological conditions of 1 year or several years. There are three common types of cycling schemes for generating the forcing using (a) observations of a single year (hereafter the single-year cycling scheme) (e.g., W. Chen et al., 2003;Lawrence et al., 2008); (b) observations of multiple years or their averages (e.g., Dankers et al., 2011;Park et al., 2013); and (c) all available observations (e.g., Barman & Jain, 2016). We refer to the last two types as the multi-year cycling scheme. Because of the large soil memories, the strategy of creating the spin-up forcing including the total spin-up length and cycling scheme, inevitably affects subsequent simulations (Schlosser et al., 2000). Therefore, systematic investigations of the impacts of the spin-up strategy on modeling regional permafrost dynamics appear indispensable, especially when dealing with regional permafrost changes in the context of global climate change.
The specific objectives of this study include: (a) exploring the efficiency and behavior of model convergence, as well as the impacts of different spin-up lengths on the subsequent simulation; (b) evaluating the impacts of cycling schemes on the resulting initial thermal fields and the consequences when simulating permafrost changes over the QTP; and (c) discussing the implications and recommend a plausible spin-up strategy for simulating regional permafrost change on the QTP.
The datasets and methodology are described in Section 2. The impacts of spin-up strategies on the resultant initial fields and simulations of changes in permafrost area are presented in Section 3. Explanations of the model convergence process, recommendations for a spin-up strategy, and implications for permafrost modeling are provided in Section 4. Section 5 concludes the paper.

Study Area and Data Sources
The QTP is bounded by 25-40°N and 75-105°E and has an area of about 2.6 million km 2 (Figure 1). It is rugged with an average elevation exceeding 4,000 m above sea level (a.s.l.). The mean annual air temperature (MAAT) is  (Tang, 2019). The distribution map is the result of the experiment SP_79-83, which forms the spin-up forcing by cycling the 1979-1983 meteorological conditions. around −2.0°C over the QTP during past decades and decreases with increasing elevation. In winter, the QTP is dominated by the cold and dry westerly, while in summer, the moisture-laden winds (the southwest monsoon) from the Indian Ocean enter into the plateau (K. Yang et al., 2014). Precipitation on the plateau varies spatially (L. Wang, Yao, et al., 2021), declining from the southeast to northwest. Annual precipitation is more than 1,000 mm in the southeast, while it is less than 50 mm in the northwest. Most of precipitation falls in the form of rain in the warm season. The average snow-depth during winter (December-February) is less than 10 cm in most locations except in the southeast high mountains (Che et al., 2008). Vegetation types change from forests in the southeast, to meadows and steppes in the center, and to deserts in the northwest, depending on temperature and precipitation conditions (Z. . The QTP is also the largest permafrost region at middle and low latitudes (G. D. Cheng & Wu, 2007). Almost half of the plateau is underlain by permafrost (Zou et al., 2017), while the rest is mainly by seasonally frozen ground. Alpine continuous permafrost is widely distributed over the Qiangtang High-plateau above 4,500 m a.s.l. In this study, we followed the permafrost continuity classification system prescribed for high-altitude regions (G. Cheng & Wang, 1982), which has the prefix "alpine" to distinguish it from the high latitude classification (van Everdingen, 1998).
Only a limited number of meteorological stations are available on the QTP and most of them are on the eastern plateau, which is mostly underlain with seasonally frozen ground (K. Yang et al., 2014). The Tanggula site (TGL) is situated on the interior plateau (91°56´ E, 33°04´ N; 5,100 m a.s.l.), with an MAAT of around −4.9°C and an annual precipitation of 400 mm in the 2000 s. This site is unique because of its comprehensive measurements of permafrost and environmental factors. It includes an active-layer monitoring system for measuring soil temperature and liquid water content at depths of 5, 10,20,40,70,105,140,175,210,245,280 and 300 cm below the ground surface, and a borehole system for monitoring ground temperature down to 34.5 m (Xiao et al., 2013). The observation began since June 2004.
Despite many available global reanalysis datasets, most of them have been found to be substantially biased on the QTP due to the presence of complex mountainous terrains and unique atmospheric circulation patterns (A. Wang & Zeng, 2012;L. Zhang et al., 2020;Zhou et al., 2021). The recently developed gridded reanalysis dataset, that is, the China Meteorological Forcing Dataset (CMFD) (He et al., 2020), is one of those widely used for the QTP studies (Han et al., 2021;T. Wang et al., 2019). It was created by fusing site-scale observations from the China Meteorological Administration with existing reanalysis data, such as Global Land Data Assimilation System (GLDAS), Tropical Rainfall Measuring Mission, and Global Energy and Water Exchanges-Surface Radiation (He et al., 2020;K. Yang et al., 2010). The dataset contains seven near-surface meteorological variables: 2 m air temperature, air pressure, specific humidity, 10 m wind speed, downward shortwave radiation, downward longwave radiation and precipitation rate. It is currently provided at 0.1° and 3-hourly for 1979-2018. According to the CMFD over the QTP, MAAT increased at a rate of 0.053°C/a from 1979 to 2010, and annual precipitation simultaneously increased from 300 to 400 to 400-500 mm in the mid-1990s ( Figure 2).
Land cover types on the QTP were retrieved from the 1: 1,000,000 Vegetation Atlas of China (Hou, 2001). Soil texture information was obtained from a multi-layered soil texture dataset that incorporates China soil properties dataset and the surveyed deep boreholes on the QTP (Wu & Nan, 2016).

Noah LSM and Modifications
The Noah LSM simulates one-dimensional water and heat exchange and transfer processes in the atmosphere-canopy/snowpack-soil continuum (F. Chen et al., 1997). Due to its concise model structure, robust parametrization, and high computational efficiency, the Noah LSM has been employed as the land scheme in the Weather Research and Forecasting Model (Ek et al., 2003) and GLDAS (Rodell et al., 2004), and has also been applied to data-scare areas, such as the QTP (van der Velde et al., 2009;Zheng et al., 2014).
In the Noah LSM, a single layer parameterization that considers fractional snow coverage and snow density change due to compaction is employed to represent the snowpack (Koren et al., 1999). The sensible heat flux is calculated based on Monin-Obukhov similarity theory, which is related to atmospheric stability (F. Chen et al., 1997). Potential evaporation is estimated by the Penman approach and constrained by the Jarvis-type resistance parameterization for actual evaporation (F. Chen et al., 1996). The surface condition is determined by the predominant land cover type, without consideration of subgrid variability.
The Noah LSM models the heat transfer in the soil following a heat diffusion equation: where, T is soil temperature; ρ, water density, and; L, latent heat. t and z represent temporal and depth dimensions, respectively. The last term of right-hand side represents heat sink/source caused by phase change. Both soil heat capacity, C s , and thermal conductivity, κ h , are functions of soil type and soil moisture (i.e., liquid water plus ice) content, θ. Soil thermal conductivity is parameterized by a normalized scheme proposed by Johansen (1975). Soil heat capacity is the arithmetic mean of the heat capacity of each soil composite weighted by corresponding volumetric content, which can be estimated as: where, the subscripts (i.e., liq, ice, air, solid) represent the liquid water, ice, air, and solid particle, respectively.
The migration of liquid water through the soil is based on Richards' equation: where, D(θ) is the soil hydraulic diffusivity, and K(θ) is soil hydraulic conductivity. Both of them are parameterized by the scheme proposed by Campbell (1974), in which the soil-relevant parameters are determined by the pedotransfer function. S(θ) represents the sink/source term to account for infiltration, evapotranspiration and phase change. Both soil heat transfer and water migration equations are solved by the numerically stable Crank-Nicolson method (Kalnay & Kanamitsu, 1988).
In the Noah LSM setting, soil water and heat are linked by the phase change and the effect of soil moisture on soil thermal properties. The convective and advective heat fluxes associated with water flow are simply neglected. The phase change process is implemented by the Clausius-Clapeyron equation, where the vapor pressure over pure ice is assumed to be in equilibrium with the soil water potential (Koren et al., 1999). Using this method, unfrozen water can be present in frozen soils, as observed in field experiments. The surface runoff scheme is based on the Simple Water Balance Model that considers spatial variabilities of infiltration capacity .
In recent years, we have improved the physical processes of the Noah LSM from three aspects in order to better represent permafrost dynamics on the QTP (H. Chen et al., 2015;Wu et al., 2018). A new thermodynamic roughness algorithm (K. Yang et al., 2008) was introduced to parameterize well the turbulent flux transfer over widespread sparse vegetation cover. The thermal conductivity parameterization scheme (Côté & Konrad, 2005) was adapted to the coarse-grained soils prevalent on the QTP. In wake of the impedance of ground ice to liquid water transport (Shoop & Bigl, 1997), we constrained the soil hydraulic conductivity by imposing an exponential empirical factor related to ice content (X. Zhang et al., 2007).
Apart from these, to account for the "heat reservoir" effect, the simulation depth was extended to 15.2 m, exceeding the average DZAA, which is about 8-10 m on the QTP (Xie et al., 2015), with a stratification of 18 soil layers to enhance the simulation of the freeze-thaw front movement . Accordingly, for numerical stability, a fixed temperature boundary was set as the lower boundary of the heat transfer equation at 40 m ( Figure 3). It was determined for each cell by using a relationship between the ground temperature at 40 m and the elevation of that cell that was derived from deep borehole observations (H. Chen et al., 2015). A free drainage boundary was used for solving the water migration equation.
The modified model has been well verified using observations at the TGL site, ground temperatures monitoring in boreholes along the Qinghai-Tibet highway, and regional permafrost maps in some surveyed areas . More details about the Noah LSM improvements and validations can be found in H. Chen et al. (2015) and Wu et al. (2018). Although there is a lack of sophisticated parameterization of snowpack and consideration of organic soils, the modifications we have made are applicable to most regions of the QTP due to the overall thin and patchy snowpack (Che et al., 2008) and thin organic soils (D. Wang, Wu, et al., 2021). According to Jin et al. (2008), the insulating effect of snowpack on the QTP is effective only when the snowpack exceeds a thickness of 20 cm on the QTP. These augmentations could effectively mitigate the cold biases in soil temperature and the overestimation of permafrost area that occurred in the original version.

Experimental Design
To test the impact of the spin-up strategy on permafrost simulation, we purposely designed two groups of experiments, one group to test varying total spin-up lengths and the other for the different cycling schemes. The same modified Noah LSM and model settings are used in all experiments, except for the spin-up forcing data which vary in each experiment. The model runs for the gridded QTP with a spatial and temporal resolution of 0.1° and 3 hr in consistency with the CMFD.
Due to lack of reliable information on soil thermal and moisture states, we have to uniformly initialize the model state variables in all grid cells using the observations at the TGL site prior to starting the spin-up runs. More specifically, the initial soil temperatures for all simulated soil layers defined at the midpoints of those layers were

of 23
linearly interpolated from the soil temperatures observed on 1 January 2008 at TGL. The volumetric soil ice contents were initialized equal to field capacity, and all other state variables used their default values. Next, the spin-up forcing was constructed using different strategies and the spin-up simulations were run. Once the spin-up is complete, the resulting state variables provide initial fields to the subsequent change simulation.
The first group was designed to investigate the impact of total spin-up length ( Table 1). The identical forcing of the first year (i.e., 1979) was cycled for varying lengths to construct the forcing data for spin-up. After spin-up, the simulation was run for the period of 1979-2010 (experimental period). In this group, we tested various lengths of spin-up, that is, 30, 100, 200, 500, 1,000 and 2,000 years, in the hope of investigating the convergence behavior of the model and determining an optimal length as a tradeoff between computational cost and accuracy.
Optimal spin-up length is also likely related to the very initial conditions provided for the spin-up simulation. We tested three variants of initialization at two locations: a permafrost location (TGL) and a seasonal frost location (Maqu), with the same temperature profile on 1 January 2008 at the TGL site and different soil moisture configurations. In the first configuration (DRY), the soil column dries out completely; in the second (LIQ), the soil layers are filled with a field capacity pore volume of liquid water; and in the third (ICE), the fraction of liquid water in LIQ turns to ice. The field capacity was predetermined by the soil texture of each soil layer in each grid cell. The test results are shown later, and based on the results, ICE was selected for the subsequent spin-up experiments.
The second group targets to investigate the impacts of cycling strategies of using a single year or multiple years to construct the spin-up forcing ( Table 2). The total spin-up lengths for the experiments were set as 2,000 years. For the multi-year cycling schemes, the total spin-up lengths were expanded accordingly to accommodate full cycles.  The mean annual air temperature (MAAT) anomalies are calculated relative to the mean of 1979-1983, which is −2.31°C. To investigate the impacts of annual atmospheric variations introduced by the single-year cycling schemes, we used the first 5 years of the available data, which are highly fluctuating, to constitute single-year cycling experiments, respectively (i.e., SP_1979, SP_1980, SP_1981, SP_1982, and SP_1983). Their experimental periods were set to 1979-2010. Compared to single-year cycling schemes, multi-year cycling schemes are theoretically better at reducing the possibility of using an abnormal year as a cycling year and are thus widely adopted in spin-up practices (Schlosser et al., 2000). We used the first five years (SP_79-83) and 10 years (SP_79-88) of forcing data as cycling years to create the multi-year cycling experiments. The experimental periods start next to the end of the cycling years ( Table 2). As a reference for comparison, SP_AVE was calculated as the average of the results of five single-year cycling experiments. Note that SP_AVE is not a real spin-up simulation, nor does it represent the truth. Rather, SP_AVE served as a reference to which the anomalies of the results from the second group of experiments were computed and analyzed. By this means, many methodological errors common in the experiments can be eliminated when working with the anomalies.
Naturally, the spin-up period should avoid to overlap with the experimental period. However, when the data length is limited, some studies even used all available data for both spin-up and simulation to maximize the experimental period (e.g., Barman & Jain, 2016). To investigate the potential impacts of overlapping spin-up and experimental simulation, we designed four experiments of multi-year cycles, including SP_79-83, SP_79-88, SP_79-98 and SP_79-08 cycling with 5, 10, 20 and 30 years of atmospheric data, respectively. The experimental periods were fixed to 1979-2010 (Table 2), which overlap with the spin-up period to varying degrees. In addition, to contrast with SP_79-08, which has an obvious warming trend during its cycling years, a detrended experiment (SP_79-08-M) was constructed using the averaged climatic conditions during this period as the cycling element.

Identification of Model Convergence
The model outputs were set at a 3 hr temporal step and were aggregated to monthly values before applying the evaluation metrics below. Referring to previous studies, for example, Z.-L. Yang et al. (1995) and Cosgrove et al. (2003), we define model convergence as the differences of monthly mean state variables that are within a threshold from year to year: where, x indicates the monthly value of a state variable at the i-th soil layer; m is month; n and n-1 represent the n-th year and the previous year, respectively. In this study, we focused on soil temperature and total soil moisture content. We tested 1, 0.1 and 0.01°C as thresholds for soil temperature, while for soil moisture content we tested 0.1, 0.01 and 0.001 m 3 m −3 The middle quantities are predefined in accordance with the accuracy of the observation system.

Quantification of the Differences in the Resultant Initial Fields
The resulting initial fields are determined as the values of state variables of the last month (December) at the end of the spin-up simulation. For the first group, we compared the resulting initial fields of soil temperature with those of SP_2000, which were found to be fully converged. Please refer to the results section which shows that.
We regarded a perfect agreement if the differences were within the accuracy range of the observation system, that is, ±0.1°C, and visually ranked the differences accordingly. For the second group, we subtracted the averaged values (SP_AVE) from the resulting initial fields of all experiments and then compared their anomalies. In visualizing the results, we also followed the same ±0.1°C divisions.
Four soil layers centered at representative depths, that is, TOP ( Mean absolute error (MAE) was employed to measure the differences of the resulting initial field with respect to the referenced field at each specific soil layer: where x and y are the values of a certain state variable at the i-th grid cell of the initial fields resulting from the reference experiment and the other experiments, and is the total number of grid cells and amounts to 25,382 in this study.

Determination of Permafrost Area
In permafrost change studies, we are often concerned with the change of permafrost area in response to climate change Obu et al., 2019). Permafrost area is calculated as the total area of permafrost cells. A cell is categorized as a permafrost type if any soil layer in this cell has maintained a maximum temperature of 0°C or below for at least two consecutive years. Otherwise, a cell is classified as a type of seasonally frozen ground if any layer in this cell has a minimum temperature equal to or below 0°C. The rest of cells are of unfrozen ground type. According to this definition, the cell type in a given year depends on the ground thermal regimes in that year and the preceding year. In view of modeling accuracy, the classification was based on aggregated monthly ground temperate to avoid potential misclassification, especially in warm permafrost regions.
In each experiment, after the spin-up simulation was completed, we continued to simulate for the experimental period driven by the real climate forcing. We compared the impact in terms of permafrost area change and degradation rate. The degradation rate is simply defined as the slope of the linear regression of simulated permafrost areas through the experimental period. For those experiments where spin-up and the experimental period did not overlap, the common phase (1990-2010) was chosen for comparison. Otherwise, the period of 1980-2010 was chosen. In addition, changes in permafrost area for the last two decades of the spin-up period were also examined. Figure 4 shows the increase in percentage of convergent cells along with total spin-up length. Overall, similar behavior is observed for both soil temperature and moisture content, regardless of the given thresholds. However, the convergence progress varies, depending on the state variable of interest as well as the given thresholds. With the less strict thresholds (i.e., 1 and 0.1°C for soil temperature, and 0.1 and 0.01 m 3 m −3 for total soil moisture), most grid cells reach convergence rapidly in the first several decades of spin-up period. The rates slow down in the next stage before convergence is reached in almost 100% cells. Given a threshold of 1°C in soil temperature and 0.1 m 3 m −3 in total soil moisture, convergence can be reached in all cells in 30 years, while for a threshold of 0.1°C in soil temperature and 0.01 m 3 m −3 in total soil moisture, in about 500 years. Given the strictest thresholds (i.e., 0.01°C for soil temperature and 0.001 m 3 m −3 for total soil moisture), after 2,000 years of spin-up simulation, the model still fails to converge in about 10% of cells. Considering the precisions of the observation instruments, the thresholds of 0.1°C for soil temperature and 0.01 m 3 m −3 for total soil moisture satisfy our study. Therefore, the experiment SP_2000 with a length of 2,000 years can be safely regarded as the reference experiment for the first group of experiments. From Figure 4, it is also evident that convergence in terms of soil moisture is faster than that of soil temperature. We will focus on convergence in terms of soil temperature in the remaining contents.

Spin-Up Length Impacts on Model Convergence
We subtracted the resulting initial thermal field of the reference spin-up experiment (SP_2000) from the fields of the first group experiments and the spatial discrepancies are shown in Figure 5. The grey areas (within ±0.1°C in discrepancy) indicate very good agreement at those locations between the benchmark experiment (SP_2000) and the testing experiments. The warm colors (orange and red) indicate positive differences whereas the cold colors (green and blue), negative differences. It is obvious that extreme discrepancies (>1.0°C or <−1.0°C as in red or blue) are rare at any depths in the study area and the discrepancies in vast areas fall within the range of ±1.0°C (in orange and green). The areas with discrepancies exceeding ±0.1°C primarily concentrate on the Qiangtang High-plateau, western and northeastern parts of the plateau, with similar spatial patterns across the experiments albeit in contrasting extent. Note that, according to Figure 4, almost all cells have reached convergence in about 500 years but in some locations, especially in the MID layer, distinct discrepancies still exist between SP_500 and SP_2000 ( Figure 5). This is due to drift in finding numerical solution . The converged values may still subject to minor changes as the computation continues.
As the total spin-up length increases (from the top down in Figure 5), the number of grid cells with detectable discrepancies (exceeding ±0.1°C, non-grey areas) gradually decreases at all depths. The colored areas have disappeared in a direction from southeast to northwest along with increasing spin-up length. In the TOP layer (the first column in Figure 5), the cells with detectable discrepancies are widespread on the northern and western parts of the plateau after a 30-year spin-up simulation. Most of them have disappeared after a 1000-year spin-up with only sporadic appearance left somewhere in the central and eastern QTP. The trends are similar for all depths as clearly presented in Table 3. Meanwhile, as a regional average measure, the MAE reduces as the spin-up length increases. Averaging the four layers, the MAE drops considerably from 0.26°C with a spin-up length of 30 years to 0.05°C with 1,000 years.
In regions with extensive occurrence of seasonal frost, such as the southern and southeastern QTP, discrepancies are rapidly diminished to an undetectable level (within ±0.1°C) even with a 30-year spin-up simulation. It corresponds to the rapid increase in the percentage of convergent grid cells at the beginning of the spin-up (Figure 4). However, in the low-elevation Qaidam Basin on the northeastern QTP, large discrepancies remain after a 30-year spin-up although seasonal frost dominates this area. It is likely connected to the fact that the initial conditions provided to the spin-up (which comes from the TGL site that is a permafrost site) are highly distinct from the real thermal and hydrological conditions in the Qaidam Basin. Nevertheless, those discrepancies are largely reduced after a 100-year or longer simulation. By contrast, areas with detectable discrepancy persist over the Qiangtang High-plateau, western and northeastern parts of the plateau even after a considerable length of spin-up simulation. These regions are predominantly underlain by permafrost as indicated by Figure 1c. This indicates that the presence of the underlying permafrost prolongs the numerical convergence and requires much longer spin-up simulation duration than at non-permafrost locations. Convergence also differs along the soil depths. Using SP_200 as an example, the shallowest layer (TOP) has the fewest areas with detectable discrepancy among the four layers, while the active layer bottom layer (MID) has the most. The 10 M and BOT layers fall in between, with similar spatial patterns of discrepancy. The positively biased areas (orange) of the interior QTP in the TOP layer transition to negatively biased areas (green) in the MID layer. The MID layer is located at a depth slightly below the permafrost table on the QTP. It hints that the existence of perennial ice in the ground (MID) may alter the spin-up behavior in permafrost regions. Compared to the shallowest soil layer, the depths within the permafrost body are more challenging for convergence and the depth of permafrost table appears to be the most challenging. We will explain the different convergence behavior in response to the same external climatic forcing in frozen soil regions later in the discussion section.
In the initialization experiments with varied soil moisture configurations at two specific sites (Figure 6), the ICE configuration (initialized with ice) appears to be superior to the other two configurations for the spin-up simulation in terms of both soil temperature and moisture at the permafrost location (TGL). Compared to TGL_DRY (assuming complete dryness), the time to reach equilibrium is reduced by about 100 years for both soil temperature and soil moisture. At the seasonal frost site (Maqu), Maqu_LIQ (initialized with liquid water) shortens the time to equilibrium the most, but the advantage of Maqu_LIQ over Maqu_DRY and Maqu_ICE is minor. All of these configurations lead to a fast reaching of equilibrium at Maqu. This finding justifies the choice of the ICE configuration to initialize the spin-up simulations in the following analysis.  Note. The mean absolute error (MAE) is calculated by regionally averaging the spatial discrepancies between the initial thermal fields resulting from the test experiments and that from SP_2000 as a reference. TOP, MID, 10 M and BOT represent different depths. The numbers in parentheses are the percentages of grid cells with discrepancies beyond the ±0.1°C range. The distributions of the discrepancies are shown in Figure 5.

Impacts of Cycling Years Scheme on the Resulting Initial Field
Unlike the first group, there is no "true" experiment in the second group. Instead, we used as reference the SP_ AVE, which is the average of the results of five single-year cycling experiments. The anomalies in the initial thermal field resulting from all experiments after a total spin-up length of around 2,000 years presented in Figure 7 are relative to the average field of SP_AVE. In general, the anomalies are pronounced and pervasive throughout the QTP and their spatial patterns appear to be remarkably distinguishable from each other at each depth. This indicates the choice of cycling years has a broad impact to the resulting initial thermal field both horizontally and vertically. The widespread presence of anomalies, as shown in Figure 7, appear to be independent of the frozen ground type.
In the near-surface layer (TOP), anomalies greater than ±1.0°C (red and blue) are found on many areas in all experiments except SP_1983 and SP_79-83. The anomalies diminish in deeper layers by turning orange and green in most experiments. It implies the meteorological conditions in the spin-up period, determined by the cycling strategy, highly affect the resulting initial field. Although many factors such as snowfall and snow cover, and rainfall in summer all influence the energy budget, air temperature generally plays the most important role in affecting thermal regimes in permafrost regions (G. . As given in Table 2  . Three initialization variants configured with different soil moisture conditions were tested. DRY: assumed completely dry. LIQ: assumed to be filled with liquid water to field capacity. ICE: assumed to be filled with ice to field capacity. with a mean higher than SP_AVE. The anomalies corresponding to the three experiments are all positive as expected. The strongest impact comes from SP_79-08, leading to the widespread of red areas (exceeding 1.0°C in anomaly) even at the lowest depth (Figure 7). It aligns with the highest MAAT in the spin-up period amidst the three, 0.80°C higher than the 1979-1983 mean. While SP_79-08-M and SP_79-08 have an identical MAAT in the spin-up period, the former has eliminated the inter-annual trend from the cycling years, resulting in much less extreme anomalies distributed in space and depth compared to SP_79-08. It thus highlights the importance of choosing representative year(s) to construct the spin-up forcing by cycling them. However, when a singleyear cycling scheme is used, it is often hard to determine a proper year given the existence of strong interannual climatic variability and usually no knowledge about the climate conditions at the starting point of simulation. Vertically, the impacts are attenuated along the depth, as demonstrated in all experiments (Figure 7). This is physically in accordance with Stefan's equation, according to which the amplitude of the soil temperature dampens exponentially with increasing depth. The regional averaged anomalies in terms of MAE apparently indicate this trend (Table 4). The anomalies are very high in the near-surface layer (TOP) and decrease sharply down to the depth of permafrost table (MID), from 1.19°C as an averaged MAE of all experiments at TOP to 0.58°C at MID. The reduction extends downward, but at a slower rate (Table 4). In terms of spatial pattern, the anomalies at TOP are distinctly different from those at MID. For example, in SP_1981, negative anomalies (blue and green) dominate at TOP whereas positive (orange), at MID. In contrast, changes in the spatial pattern from MID downward are not as evident.
It is interesting when we look at SP_79-83 and SP_1983 results shown in Figure 7, the former ending with the year of 1983 and the latter cycling over the same year. The two have a highly similar spatial pattern of anomaly at TOP. Such similarities are gradually weakened at the lower depths. In SP_1983 the widespread negative anomalies persist through the depth, but in SP_79-83 most negative anomalies at TOP have disappeared at 10 M and BOT. As measured by MAE, the anomalies in SP_1983 are almost twice those in SP_79-83 in all layers, except for the TOP layer where they are close between the two experiments. It well demonstrates the variable responses at different depths to external climatic forcing. The thermal fields of shallow soils are susceptible to short-term meteorological conditions (e.g., 1983) and those of deeper soils tend to mainly respond to long-term climatic variations. Since SP_AVE provides a field averaged from the five single-year cycling experiments and represents an average state, it is not surprising that it matches well with SP_79-83 in the deep layers (10 M and BOT). Therefore, it makes the multi-year cycling scheme preferable to the single-year cycling scheme, as the interannual variability is likely more representative for the starting point of simulation than a single year, while short-term effects on the shallow depths resulting from the cycling years can be easily eliminated in a follow-up simulation where actual climate data drive the model.

Spin-Up Length Impacts on the Simulation of Permafrost Area
A spin-up simulation provides appropriate initial states of the model variables to subsequent simulation. Thus, the use of different spin-up strategies indirectly affects the results of change simulations. Figure 8 presents the change in permafrost areas over time in response to various spin-up strategies. Figures 8a and 8b show changes in plateau permafrost area over the last two decades of the spin-up period and Figures 8c and 8d changes over the follow-up experimental period.
In the first group of experiments on varying spin-up lengths, the change pattern of SP_30 differs profoundly from the others (Figure 8a). At the beginning (20 years before the end of the spin-up period), SP_30 has simulated many more permafrost areas before rapidly converging to a stable value, close to that simulated by other experiments. This change trend demonstrates well a convergence process during a spin-up simulation that gradually narrows the discrepancy, as the TGL soil temperature profile has been uniformly applied to all cells on the QTP as the initial states, which probably deviates from the real thermal field. The other experiments overlap with each other regarding the permafrost area in the last 20 years until the end of spin-up (Figure 8a). Thus, if we treat permafrost area as a convergence criterion, we are likely to be misled as if a 100-year length of spin-up simulation be sufficient for convergence. However, as shown in Figure 5 and Table 3, the final thermal regimes on the QTP, especially in permafrost areas, produced in SP_100 are far from stabilized. This finding indicates that, although changes in permafrost area are of the greatest concern, permafrost area is simply inappropriate for use as a criterion for determining model convergence status.
Despite no obvious areal differences detected in the experiments on total spin-up lengths ≥100 yeas, the unperceived discrepancies in the resulting initial fields ( Figure 5 and Table 3) come to exert influence on the follow-up simulation informed by actual climatic conditions (Figure 8c). In the first decade (1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989) (Figure 8c). However, the impacts emerge and continue to increase as the simulation progresses. Compared with SP_2000, SP_30 with the shortest spin-up length severely underestimates the permafrost area at the end of the simulation (2010), with a discrepancy up to 115.2 × 10 3 km 2 , which approximates 11% of the area simulated by SP_2000. Meanwhile, SP_500 and SP_1000 show good agreement with SP_2000. Therefore, sufficient spin-up lengths are critical to simulate changes in permafrost area and inadequate one can lead to a large underestimation of permafrost area.
Based on those simulations, the calculated permafrost degradation rates are shown in Table 5. It is evident that the degradation rate decreases with growing spin-up length. In other words, insufficient spin-up lengths tend to overestimate the permafrost degradation rate. SP_30 can bias the rate by as much as 3.98 × 10 3 km 2 /a, approaching 86.3% of the SP_2000 rate. The rates from SP_500 and SP_1000 are close to the SP_2000 rate with discrepancies of less than 10%. Considering its overall performance in the resulting initial field and the follow-up simulation, SP_500 provides an acceptable accuracy in estimating permafrost area change while offering computational advantages over SP_1000 and SP_2000.

Impacts of Cycling Years Scheme on the Simulation of Permafrost Area
The impacts of cycling schemes on the simulated permafrost area are pronounced, and the changes over time appear to be more complex than those due to varying spin-up lengths. In all experiments on single-year cycling  Note. The period of 1990-2010 is used for analysis to ensure no overlap between the cycling years and the experimental period in SP_79-83 and SP_79-88.

Table 5
Simulated Rates of Permafrost Area Degradation schemes, constant permafrost areas are estimated in the last two decades before the end of a 2000-year spin-up simulation (Figure 8b). By contrast, the experiments on multi-year cycling schemes with a similar total spin-up deliver oscillating results of permafrost area. Among them, the experiments cycling with the first 5-year (SP_79-83) and 10-year (SP_79-88) of the forcing data produce dynamically-stable permafrost areas with a fixed period. For SP_79-98 and SP_79-08, Figure 8b shows only part of the intra-period variations of the last cycle because the 20-year display window falls within a single cycle (20 years, SP_79-98; 30 years, SP_79-08).
Assuming SP_AVE as the baseline, the positions of the cycling experiments that lie above or below SP_AVE (Figure 8b) are roughly opposite in sign to the MAAT anomalies of the cycling years with respect to the 1979-1983 mean (Table 2). This contrasting pattern persists through the subsequent experimental period and is still visible through 2010 (Figure 8d), suggesting that the effects of cycling schemes are sustained for at least decades after the spin-up simulation ends. Using the initial fields provided by these spin-up simulations, the follow-up simulations show distinctly different patterns of area change through the experiment period and different magnitudes of permafrost areas at the end of the period (2010), as shown in Figure 8d. The largest discrepancy in simulated permafrost area in 2010 is 173.5 × 10 3 km 2 , about 16% of the area in the SP_AVE. Compared to the impacts caused by different spin-up lengths, the application of different cycling schemes in the spin-up simulation leads to more profound impacts on the simulation of permafrost area.
In terms of estimated permafrost degradation rate, the discrepancy caused by different cycling schemes reaches up to 3.26 × 10 3 km 2 /a during 1990-2010, which was specially chosen to ensure no overlap between the cycling years and the experimental periods (Table 5). The difference is 2.75 × 10 3 km 2 /a among five single-year cycling experiments, while the value is only 0.53 × 10 3 km 2 /a for the multi-year experiments (i.e., . This is a good evidence of the robustness of multi-year cycling schemes, which can reduce the large uncertainty of single-year scheme caused by strong atmospheric variations. In the case of allowing overlap between cycling years and experimental periods in the simulations, the largest difference in the rates of permafrost area degradation during 1980-2010 reaches 3.37 × 10 3 km 2 /a ( Table 5). As for those single-year cycling schemes, the difference can be up to 2.87 × 10 3 km 2 /a, while it is 2.66 × 10 3 km 2 /a for those multi-year cycling schemes. SP_79-83 and SP_79-88 indicate similar rates with a small difference of 0.65 × 10 3 km 2 /a. As for SP_79-98 and SP_79-08, they have the smallest degradation rates, due to a large increase in simulated permafrost area in 1980-1985 that offsets the declining trend after 1990 (Figure 8d). The large increase in the early stage of the simulation is believed to be closely connected to the apparent warming trend introduced in the long cycling years in both schemes. We will discuss this issue in a later section. In comparison, SP_79-08-M is a detrended version of SP_79-08 with interannual variability removed from the cycling years of SP_79-08. However, due to overall higher MAAT in the cycling years of SP_79-08-M, the sharp increase in simulated permafrost area is still present in the early stage of the simulation (Figure 8).

Model Convergence Behavior in Cold Regions
In permafrost regions, model convergence rates vary with depth ( Figure 5 and Table 3). In theory, shallow soils respond quickly to atmospheric variations, so model converges relatively more quickly for soil temperature in the TOP layer. At depths, such as the 10 M and BOT layers, despite the suppressed thermal responses to external forcing, the variations in soil temperature are minor in space in comparison with those in the shallow soils, making the TGL profile well approximating to their equilibrium states. As a combined result, relatively rapid convergence can still be achieved in many permafrost cells even at greater depths. Meanwhile, similar spatial patterns of discrepancy are found between the layers of 10 M and BOT as a result of the small soil temperature gradient between the deep soils.
The MID layer is generally near the bottom of the active layer in permafrost regions on the QTP, where phase changes substantially alter the partitioning of energy and water, while they do not respond effectively to atmospheric variations, unlike shallow soils. More importantly, drastically weakened permeability of ice-rich soils near the permafrost table impedes soil moisture migration . This eventually leads to complicated thermal and hydrological interactions there and a prolonged convergence of soil temperature at the freezing and thawing front that moves slowly downwards to the bottom of the active layer. Such a disequilibrium in turn affects the thermal states of adjacent soil layers through conduction due to thermal gradients, delaying convergence in these layers. In non-permafrost regions on the QTP, the MID layer is located in thawed soils. The maximum freezing depth in regions of seasonally frozen ground is less than 2 m in fine-grained Quaternary deposits (Luo et al., 2020), which is less than the maximum thawing depth in permafrost regions, which is averaged at 3 m on the QTP . It also helps explain the difficulty in model convergence in permafrost regions compared to that in regions of seasonally frozen ground and unfrozen ground.
In addition, some previous studies emphasized that the spin-up simulation should be initialized with snow cover free or thawed conditions in order to save the overall length of spin-up (e.g., Mitchell et al., 2005;Shrestha & Houser, 2010). However, as the results of our study show, these treatments may not play a vital role in spatial modeling because the freeze-thaw processes occurring at the bottom of the active layer dominate the convergence of the model simulation. On the other hand, appropriate initialization for spin-up simulation helps reduce the total length of sufficient spin-up ( Figure 6); that is consistent with previous studies (e.g., Shrestha & Houser, 2010;Z.-L. Yang et al., 1995).
As shown in Figure 8c, spin-up simulations with inadequate total spin-up lengths can substantially bias the estimates of permafrost area in the late stages of follow-up simulations forced by real climate data. This is unusual since the impacts of the initial fields resulting from the spin-up simulations should fade out as the simulation proceeds. It is partially related to the disequilibrium in the distribution of ground ice by the end of the insufficient spin-up simulations. Small perturbations in the initial fields do not immediately transform the frozen ground type, even under the real climate forcing. However, after the first decade of warming as shown in the real forcing (Figure 2), the perturbations begin to take effect in the transition areas on the QTP with warm permafrost.
In recent decades, permafrost on the QTP has degraded, mostly in the areas with high temperature permafrost and peripheral to the alpine continuous permafrost region, becoming seasonally frozen ground (Figure 9a). In order to validate the permafrost distribution simulation, we also included a map from a previous study (i.e., Zou et al. (2017), Figure 9b), which was compiled using a remote sensing approach and shows high similarity to the SP_79-83 map in terms of spatial distribution patterns. In our experiments, many cells in disequilibrium state are consistently located in the transition areas ( Figure 5). Given their extremely unstable thermal states, type conversion is more likely to occur even if their thermal states are affected by small perturbations.

The Longer the Better?
Apart from the computational costs, the longer the spin-up simulation runs, the more stable the model is. This is intrinsic to the nature of an LSM based on the equations for heat and water transport. However, a spin-up simulation over 2,000 years is computationally expensive. With a threshold of 0.1°C for soil temperature and 0.01 m 3 m −3 for total soil moisture content, convergence in all cells can be achieved in 500 years. In terms of permafrost area, a 1000-year spin-up simulation produces results closer to those of the 2000-year spin-up than a 500-year spin-up, but with minor gains toward the 2000-year results. We thus recommend that a 500-year spin-up would suffice for a good compromise between the discrepancies and the computational burden.
This length could be specific to the modified Noah LSM and the QTP. For other LSMs, experiments similar to this study are still needed for optimizing the spin-up length. Furthermore, as can be seen in Figure 8c, the discrepancies increase as the simulation proceeds to the late phase. It provides a clue that for change simulations conducted over a very long period of time, for example, over 100 years, it will be more suitable to use a longer spin-up length. It is especially true when an LSM is applied to simulating future permafrost change over a period from the present to the end of this century.
On the other hand, when long cycling years are used to create the spin-up forcing, there is a risk that it will be less representative of the climate conditions at the beginning of the simulation and also introduce a climate change trend to the spin-up simulation. As a result, there may be considerable discrepancies in simulation results. On the QTP, the 1979-1998 (corresponding to SP_79-98) and 1979-2008 (SP_79-08) periods are warmer than the 1979-1983 mean by 0.42 and 0.80°C, respectively, and even warmer than any single year in the first five years (1979)(1980)(1981)(1982)(1983) ( Table 2). The large increase in permafrost area in the early simulation period (1980)(1981)(1982)(1983)(1984)(1985) in the two experiments with long cycling years contradicts the decrease at the later stages, as well as the results from other experiments (Figure 8d). Therefore, despite the advantages over the single-year scheme, it is critically important to avoid use of excessively long cycle lengths in the multi-year cycling scheme, especially when the cycling years have to overlap with the simulation period.
Some spin-up studies also use an average climatic condition over a period of time to constitute the spin-up forcing (e.g., Rodell et al., 2005), which is in essence a single-year cycling scheme but has the advantage of being less likely to use anomalous annual conditions. According to our analysis, it may risk losing the fluctuating nature of the climate forcing that progressively affects deep soils during hundreds of cycles in the spin-up simulation compared to a multi-year cycling scheme ( Figure 7). However, we need further independent test to determine if there is a benefit to maintaining interannual variability. In addition, an average climatic condition over a warming period leads to a higher mean temperature, so such a spin-up simulation also leads to significant biases in the resulting initial fields.
It is worthy to note that it is generally impossible to reproduce a real initial field in a permafrost region using the dynamic equilibrium spin-up method as used in this study, in particular when the thermal states of permafrost are subject to rapid changes. Permafrost is a product of long-term interactions between multiple systems and it responds with a lag at depths to rapid regional climate change. In fact, thermal states of permafrost have never been in equilibrium with concurrent climate. However, in a data-scare region, such as the QTP, it is practical to perform spin-up simulation using the dynamic equilibrium method given the available forcing dataset. Until any longer climate dataset is available for the QTP, we can roughly assume a dynamic equilibrium with a sufficient cycling with the early years in the 1980s. To this end, our study calls for the establishment of longer-span, high-resolution forcing datasets on the QTP, as has been done in western Canada by Asong et al. (2020). Therefore, we recommend a spin-up strategy of a 500-year cycling with the first 5-10 years of forcing for regional modeling of permafrost dynamics on the QTP.

Implications on Simulations of Changes in Permafrost Area
Existing simulation studies has reported different rates of permafrost area degradation on the QTP over the last decades. Guo et al. (2012) presented a plateau permafrost degradation rate of 7.60 × 10 3 km 2 /a during 1981-2010 based on the Community Land Model coupled with the RegCM3 regional circulation model. The degradation rate was later updated to be as much as 9.20 × 10 3 km 2 /a in Guo and Wang (2013), who simulated with the same LSM but ran with a high-resolution forcing data. Based on the geographically weighted regression method, Ran et al. (2018) estimated a degradation rate of up to 9.52 × 10 3 km 2 /a on the QTP. However, another work (T. Wang et al., 2019), also using an empirical method, the logit regression, indicates a much lower degradation rate of 6.60 × 10 3 km 2 /a during 1980-2010. The maximal difference between those studies is up to 2.92 × 10 3 km 2 /a.
Many factors affect the estimate of permafrost degradation rate in those studies. Our work demonstrates from a new perspective that different spin-up strategies can lead to uncertainty even when an LSM is forced by the same meteorological data and with the same parameters. The difference of permafrost degradation rate during 1980-2010 caused by insufficient spin-up length can reach up to 3.98 × 10 3 km 2 /a. The largest discrepancy caused by different cycling schemes reaches 3.37 × 10 3 km 2 /a. These differences are comparable with the uncertainty (2.92 × 10 3 km 2 /a) associated with the existing studies based on different modeling methods and data, and even approach half the magnitude of the permafrost degradation rate itself. It means for all LSM-based simulations of regional permafrost change that do not fully account for the effects of spin-up strategies, the estimated degradation rates may lose their statistical significance.
In practice, the spin-up length for a regional study is usually determined according to the guidance from modeling studies at available sites with observations. As demonstrated in this study, the lengths that support sufficient spin-up simulation vary across space ( Figure 5). Hence, it is critical to perform a sensitivity analysis on spin-up length variation prior to spatial modeling of permafrost dynamics based on an LSM.
Some previous studies also assume the impacts of cycling schemes can be eliminated in a short period (Y. Zhang et al., 2008), albeit without a prior knowledge. They have chosen to discard the results of several years at the beginning of the simulation following a spin-up simulation in hopes of further minimizing the uncertainties associated (Barman & Jain, 2016;Burke et al., 2013). Unfortunately, as shown in this study, the impacts on the permafrost area are still evident by the end of the simulation (Figure 8d). Using the same set of meteorological data, the discrepancies in shallow soils as caused by different spin-up strategies are possible to be removed in the starting several years, while the discrepancies in deep soils remain for decades as a result of the memory strengthened by permafrost. Therefore, when conducting a regional simulation study of permafrost dynamics using any LSM configured with deep soil profiles, a proper spin-up strategy with a long enough total spin-up length and an appropriate cycling scheme should be carefully determined to best curtail the uncertainties.

Conclusions
This study evaluated the impacts of different spin-up strategies configured with total spin-up length and cycling years scheme (including single-year cycling and multi-year cycling options) on the modeling of regional permafrost dynamics on the QTP by employing the modified Noah LSM. The immediate impacts on the resulting initial fields of soil temperature were examined at four representative depths, that is, near ground surface (TOP, 0.025 m), near the permafrost table (MID, 3.7 m), near the DZAA (10 M, 10.2 m), and the bottom of the soil column (BOT, 14.2 m). We also explored the long-term impacts on the subsequent simulations driven by real meteorological forcing with respect to permafrost area change and degradation rate. We reached the following conclusions: A prolonged spin-up simulation is required at locations underlain by permafrost. It is relatively more difficult to achieve convergence at the depth of permafrost table than at other depths. The impacts of total spin-up length on the resulting initial field of soil temperature are relatively small. However, the small perturbations in the warm permafrost underlying transition zones also result in significant differences in estimates of permafrost area change and degradation rate in ensuing simulations. Unlike the total spin-up length, different cycling years schemes have a considerable effect on the resulting fields of soil temperature, the effects of which are profound and persist throughout the entire follow-up simulation period.
Inadequate total spin-up length leads to an underestimation of permafrost area and an overestimation of degradation rate in the follow-up simulation over a period of several decades. The degradation rates based on experiments with varying total spin-up lengths differ from each other by as much as 3.98 × 10 3 km 2 /a. Meanwhile, the difference in degradation rate due to different cycling schemes is up to 3.37 × 10 3 km 2 /a. They are comparable with the uncertainty (2.92 × 10 3 km 2 /a) associated with the existing studies based on different modeling methods and data, and even approach half the magnitude of the degradation rate itself. Therefore, a proper spin-up strategy is essential for curtailing uncertainties in LSM applications for regional permafrost modeling and improving comparability among these studies.
For simulations of permafrost dynamics on the QTP using the modified Noah LSM, a 500 year spin-up length is an acceptable trade-off between computational cost and convergence on most grid cells. It is not recommended to choose single-year cycling schemes if the representativeness of the chosen year does not reliably reflect real climate conditions in the starting stage of the simulation. A multi-year cycling scheme is generally preferred, but an overlong length of cycling years should be avoided to prevent the introduction of climate change trends in the