## 2 Materials and Methods
### Model Description and Development
is the density of ice, The density of ice is \(\rho_{ice}\) (917 kgm\({}^{-3}\)), \(\theta_{ice}\) (m\({}^{3}\)m\({}^{-3}\)) is the partial volume of ice content, and \(L_{f}\) (0.3336 \(\times\) 10\({}^{6}\) Jkg\({}^{-1}\)) is the latent heat of fusion. (2) Vegetation dynamic module  
After registering a change in the biomass for each plant type, this module calculates the gross primary productivity (GPP), the net primary productivity (NPP), and the net ecosystem exchange (NEE) in the vegetation coverage and type. The vegetation productivity calculation includes both the total photosynthesis and the breathing per hour. Vegetation biomass is the sum of the root system biomass and the plant leaves. The GPP calculation [20] is:  
\[GPP=\int A_{g}dt \tag{5}\]  
and the NPP calculation is:  
\[NPP=(1-\eta)\int\big{(}A_{g}-R_{\text{leaf}}-R_{\text{stem}}-R_{root}\big{)}dt \tag{6}\]  
where \(A_{g}\) (CO\({}_{2}\)/m\({}^{2}\cdot\)s) is the net photosynthetic rate of the canopy; \(R_{leaf}\) (CO\({}_{2}\)/m\({}^{2}\cdot\)s) is the leaf respiration rate; \(R_{stem}\) (CO\({}_{2}\)/m\({}^{2}\cdot\)s) is the respiration rate; and \(R_{root}\) (CO\({}_{2}\)/m\({}^{2}\cdot\)s) is the root respiration rate. The rate of photosynthesis is simulated according to the Farquhar equations [36]. In accord with [PERSON] [37], \(\eta\) is 0.33 or 0.3, which is the proportion of carbon loss in tissue construction. (3) Vegetation Phenology Module  
In the phenological module, it is assumed that deciduous plants begin to germinate in winter when the growth of a particular plant functional type exceeds a certain threshold. This threshold is defined by the day when two conditions are met: (1) the average temperature (_Avt10_) is greater than 5 \({}^{\circ}\)C for 10 consecutive days and (2) at least 100 days have experienced temperature increases greater than 0 \({}^{\circ}\)C (_GDD0_). When _Avt10_ is less than 5 \({}^{\circ}\)C,the vegetation enters the withered yellow period [38]. This threshold day forms the basis of the leaf area index (LAI) calculations [20]:  
\[LAI=LAI_{\text{max}}\times\text{f}(T) \tag{7}\]  
\[\text{f}(T)=\left\{ \begin{array}{c}(GDD_{0}-T_{\text{b}})/b\ GDD_{0}>T_{b}\ \text{and}\ avt10>5\\ 1\\ \text{f}(T)-1/dfac\ avt10<5\end{array} \right. \tag{8}\]  
where \(T_{b}\) is the threshold of \(GDD_{0}\) in the vegetation green-back period, \(b\) is the coefficient of controlling the leaf spreading speed, and _dfac_ is the coefficient of controlling the defoliation speed. (4) Modified soil erosion parameterization scheme  
Numerous studies have pointed out that soil erosion and the erosion of poorly vegetated land are accompanied by changes in soil nutrients, structure, moisture, texture, etc. The soil microbial community will also be transformed during the restoration of soil erosion and erosion-damaged land [39; 40]. As the main driver of soil organic carbon transformation, microbial community changes in erosion and revegetation systems largely determine the rate of organic carbon mineralization and fixation and the microbial contribution to soil organic carbon. The microbial communities of erosion and revegetation systems largely determine the rate of organic carbon mineralization and fixation and the apparent functional redundancy of microorganisms in the soil. Therefore, we used RUSLE to improve the soil erosion process module. The expression of general soil loss equation is:  
\[A=R\cdot K\cdot LS\cdot C\cdot P \tag{9}\]  
where \(A\) is the soil loss (kg km\({}^{-2}\)), \(R\) is the rainfall erosivity factor (J m\({}^{-2}\)), \(K\) is the soil erodibility factor (dimensionless), \(LS\) is the slope length and gradient factor (dimensionless), \(C\) is the vegetation cover and management factor (dimensionless), and \(P\) is the soil and water conservation measures factor (dimensionless).