the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 25 Jun 2025
Ensemble numerical simulation of permafrost over the Tibetan Plateau from Flexible Permafrost Model: 1950–2023
Abstract. Permafrost remains a largely subsurface phenomenon, and understanding its dynamics as well as influences under a warming climate heavily relies on numerical simulations. However, this task presents significant challenges as the state-of-the-art land surface models are weak in their ability to represent permafrost processes. In this study, we introduce a new land surface scheme specifically designed for permafrost applications, the Flexible Permafrost Model (FPM). This model serves as an adaptable framework for implementing innovative parameterizations of permafrost-related physics. The FPM accounts for both vertical and lateral heat flow at and below the soil surface, while simultaneously resolving the land-atmosphere energy exchanges through comprehensive treatment of radiative balance and turbulent flux dynamics. We simulate the ground thermal regime and test the model with a network of permafrost measurements across the Tibetan Plateau.
Our result yields root mean square error values of 1.1 m for the thickness of the active layer and 1.5 °C for the mean annual ground temperature of permafrost. We estimate that the current extent of permafrost (2010–2023) on the Tibetan Plateau is approximately 1.15 ± 0.02 × 106 km2, which aligns closely with published estimates. Long-term simulations indicate that the permafrost temperature has increased by 0.26 °C since 1980 with a decreased area of 13.2 × 104 km2 (∼10.5 %). These ensemble simulations provide valuable information on the dynamics of permafrost over the Tibetan Plateau. Furthermore, our findings suggest that current land surface models, which utilize shallow soil columns, are insufficient for permafrost simulations over the Tibetan Plateau due to the typically deep active layer (that is, 2.88 ± 0.95 m by mean) and may not be suitable for future projections.
- Preprint
(4171 KB) - Metadata XML
-
Supplement
(6 KB) - BibTeX
- EndNote
Status: open (until 21 Sep 2025)
-
RC1: 'Comment on egusphere-2025-1828', Rui Chen, 14 Aug 2025
reply
The comment was uploaded in the form of a supplement: https://egusphere.copernicus.org/preprints/2025/egusphere-2025-1828/egusphere-2025-1828-RC1-supplement.pdf
-
RC2: 'Comment on egusphere-2025-1828', Anonymous Referee #2, 21 Aug 2025
reply
General Comment
In their study, Sun and Cao, present a new permafrost model, evaluate its performance at certain locations and apply it to the whole Tibetan plateau. They compare observations with simulations forced with local weather station records and with large scale reanalysis datasets. They discuss what makes the model perform better or worse and presents results on the evolution of Tibetan permafrost since the 80s. The model includes Surface Energy Balance calculation but its calculation and its coupling to the energy budget of the soil column is, from my understanding, either poorly described or problematic in its design (see my Important Comments). It resolves heat conduction in the ground with effective heat capacity, freezing curves but the water content of each cell is static (no infiltration or upward suction via evaporation and matrix potential). At the surface it also includes a snowpack module that does not consider snow mass balance.
The model claims to be flexible and to propose novel parameterizations, but according to me the flexibility is not explained or demonstrated (see my Important Point on L58) and I do not see novel parameterizations. Additionally, I do not see what are the new possibilities that this model offers compared to already existing models. I think the study needs to better present what makes it original and motivated its development. A detailed summary of its strength and weakness compared to other models would help understand the motivations for its development.
Also, in my view, the study might be more aligned with the scope of GMD than TC, given that its primary focus is on presenting a new geoscientific model. Typically, I would expect TC publications to emphasize scientific results or insights directly related to the cryosphere, whereas GMD is intended for studies of this nature. That said, since the manuscript has already passed the initial editorial screening and the editor has decided to proceed with the review process, I take it that its placement here is considered appropriate.
Altogether, for now, I have the feeling that the model in itself is not particularly novel/needed by the community (but I am happy to be proven wrong) and, unless I am mistaken, its description includes important flaws that needs to be addressed. Therefore as it is, I recommend major revisions to address those crucial points. I have not provided detailed comments on the rest of the manuscript at this stage, as I believe the major issues outlined above should be addressed before a more thorough evaluation is meaningful. Overall, I did not identify major issues with the model setup, validation, results, or discussion. However, I find the robustness of the results at the Tibetan scale questionable, given the discrepancies between observations and simulations when using the reanalyses.
Important comments
L30-43
When a new model is published, it becomes part of the existing landscape of models, and it is important to explain how this new model positions itself relative to the existing ones in order to justify its relevance. I think the study should be more thorough in this regard and more exhaustive regarding which permafrost models are used in Tibet to study what. The Introduction should be expanded in consequence (see my comments about L42 and L43).
L58
“The application of FPM with lateral heat is provided in Sun et al. (2023).” This is very confusing to me for two reasons. First, vertical 1D models and cross sectional 2D models are usually very different types of models that implies different formulation of their physical equations, different numerical schemes for their resolution and different types of upper boundary conditions. The present study mentions surface energy balance calculation whereas Sun et al. (2023) forced their model with ground surface temperatures. So I do not understand how a 1D model can become a full 2D model (not a coupling of 1D simulation together).
Second, the given explanation is confusing. Sun et al. (2023) says “A 2D heat conduction model developed by Ling and Zhang (2004a) was used to simulate the permafrost thermal regime.” So if it is the same model in 1D and 2D, is the present study actually presenting a new model or is it presenting improvements brought to a model published initially in 2004? Since the goal of this study is to present a new model, this kind of question need to clarified to consider publication.
L62-104
“A physically-based surface energy balance scheme for different land surface cover types with varying snow regimes and properties was coupled to FPM, and was formulated as:
(1−α)Qsi +Qli +Qle +Qh +Qe +Qc = Qm”
Major problem here. The study intends to describe the ground surface SEB (because we are in part 2.1 of a study presenting a permafrost model, with no mention to the snow scheme yet). Yet it looks like a description of a snow SEB, as evidenced by the expression “varying snow regime” and by the fact that the equation calculates energy for melt, which is inappropriate for a ground surface scheme. For a ground surface scheme, we want the SEB to give us access to the dE/dt of the surface so that it can force the heat diffusion/advection in the soil column. If we do not have that, we do not have the coupling between the climate and the temperature in the ground. This problem persist over the whole section 2.1. These are very important aspect of the model description that need to be carefully addressed so that we can understand what we are talking about.
L78 and 86
I have 2 problems with this quantification of the turbulent fluxes. First, if you take Preistley and Taylor for Q_e, then you have to take Q_h as the residual of the available energy for turbulent fluxes, which is:
Q_h=(1-α ∆/(∆+γ))×(Q_n-Q_g )
Otherwise your SEB is not energy conservative. Preistley and Taylor considers that the energy available for both turbulent fluxes is Q_n - Q_g and provide a formula to find what fraction of that available energy goes to the latent flux. So, in order to have a consistent surface energy balance budget, Q_h has to be the complementary fraction of Q_n - Q_g, not another formula based on another theory.
Second, if you use Preistley and Taylor you will assume that the sum of the terms of your SEB equals 0 and not the energy variation of the soil surface. Otherwise you’d have to work with a modified version that would look like:
Q_e=α ∆/(∆+γ) (Q_n-Q_g-(∂E_surf)/∂t)
But then you would not be able to calculate together the turbulent flux and the energy variation of the surface (one too many unknown). Yet, it is through the energy variation of the surface that you can couple the SEB and the energy budget at depth in the ground (because the energy budget of the surface will drive the one of the subsurface). Therefore, this whole SEB description is very confusing to me and absolutely need to be fixed. For now I cannot understand the coupling with the subsurface.
Specific comments
L14: “shallow soil columns” indicate typical depth
L36: Wrong reference with Lan et al. 2025 here ? It is a review of the strength and weakness of a climate reanalysis datasets, it does not present a permafrost model.
L43: I am surprised the authors do not mention other models used to study Tibetan permafrost like the GBEHM model that is greatly used within the Chinese community (Fang et al., 2025; Gao et al., 2018; Qin et al., 2017; Shi et al., 2020; Wang et al., 2023, 2018; Wang and Gao, 2022, 2025; Yang et al., 2023b, a), or the recent DHTC model (Linmao et al., 2024) and FLEXTopo-FS model (Gao et al., 2022).
Eq13: I would avoid writing the equation, it gives the impression that your are describing how the model works during simulations whereas you just used the equation once for the evaluation of static parameters. I would rather just state the values of k and C in the text and say that they were calculated based on the ref you mention.
Fig.2: “Comparison of simulated active layer soil temperature with time series at the synthesis sites.” What is the methodology here? The active layer can be pretty deep. You averaged the temperature over the whole the active layer at a daily time step? Also what do the red and blue numbers with and without parenthesis correspond to? I assume from the text that it is the bias ? It should be written in the caption for more clarity.
Fig.3: The methodology described in the caption is hard to understand, please elaborate more.
Phrasing and typos
L214: Align
L218: “via the processes of latent heat and soil moisture” the reader understands, but these are not processes, please rephrase.
L221: “The remote-sensing datasets are different in temporal coverage, so we use the climatology to represent the long-term conditions.” I think the phrasing “use the climatology” can be improved.
Fig.2: “The daily soil temperature present are averaged”, syntax problem.
References
Fang, P., Wang, T., Yang, D., Yang, J., and Tang, L.: Permafrost Degradation and Concomitant Hydrological Changes Dominated by Anthropogenic Greenhouse Gas Emissions in the Northeastern Tibetan Plateau, Geophysical Research Letters, 52, e2024GL113679, https://doi.org/10.1029/2024GL113679, 2025.
Gao, B., Yang, D., Qin, Y., Wang, Y., Li, H., Zhang, Y., and Zhang, T.: Change in frozen soils and its effect on regional hydrology, upper Heihe basin, northeastern Qinghai-Tibetan Plateau, Cryosphere, 12, 657–673, https://doi.org/10.5194/tc-12-657-2018, 2018.
Gao, H., Han, C., Chen, R., Feng, Z., Wang, K., Fenicia, F., and Savenije, H.: Frozen soil hydrological modeling for a mountainous catchment northeast of the Qinghai–Tibet Plateau, Hydrology and Earth System Sciences, 26, 4187–4208, https://doi.org/10.5194/hess-26-4187-2022, 2022.
Linmao, G., Genxu, W., Chunlin, S., Shouqin, S., Kai, L., Jinlong, L., Yang, L., Biying, Z., Jiapei, M., and Peng, H.: Development of a modular distributed hydro-thermal coupled hydrological model for cold regions, Journal of Hydrology, 644, 132099, https://doi.org/10.1016/j.jhydrol.2024.132099, 2024.
Qin, Y., Wu, T., Zhao, L., Wu, X., Li, R., Xie, C., Pang, Q., Hu, G., Qiao, Y., Zhao, G., Liu, G., Zhu, X., and Hao, J.: Numerical Modeling of the Active Layer Thickness and Permafrost Thermal State Across Qinghai‐Tibetan Plateau, JGR Atmospheres, 122, https://doi.org/10.1002/2017JD026858, 2017.
Shi, R., Yang, H., and Yang, D.: Spatiotemporal variations in frozen ground and their impacts on hydrological components in the source region of the Yangtze River, Journal of Hydrology, 590, 125237, https://doi.org/10.1016/j.jhydrol.2020.125237, 2020.
Wang, T., Yang, D., Yang, Y., Zheng, G., Jin, H., Li, X., Yao, T., and Cheng, G.: Pervasive Permafrost Thaw Exacerbates Future Risk of Water Shortage Across the Tibetan Plateau, Earth’s Future, 11, e2022EF003463, https://doi.org/10.1029/2022EF003463, 2023.
Wang, X. and Gao, B.: Frozen soil change and its impact on hydrological processes in the Qinghai Lake Basin, the Qinghai-Tibetan Plateau, China, Journal of Hydrology: Regional Studies, 39, 100993, https://doi.org/10.1016/j.ejrh.2022.100993, 2022.
Wang, X. and Gao, B.: Future hydrological changes of the upstream basin and their effects on the water level of the Qinghai Lake, Journal of Hydrology: Regional Studies, 59, 102425, https://doi.org/10.1016/j.ejrh.2025.102425, 2025.
Wang, Y., Yang, H., Gao, B., Wang, T., Qin, Y., and Yang, D.: Frozen ground degradation may reduce future runoff in the headwaters of an inland river on the northeastern Tibetan Plateau, Journal of Hydrology, 564, 1153–1164, https://doi.org/10.1016/j.jhydrol.2018.07.078, 2018.
Yang, J., Wang, T., and Yang, D.: Divergent responses of permafrost degradation to precipitation increases at different seasons on the eastern Qinghai–Tibet Plateau based on modeling approach, Environ. Res. Lett., 18, 094038, https://doi.org/10.1088/1748-9326/acf05c, 2023a.
Yang, J., Wang, T., Yang, D., and Yang, Y.: Insights into runoff changes in the source region of Yellow River under frozen ground degradation, Journal of Hydrology, 617, 128892, https://doi.org/10.1016/j.jhydrol.2022.128892, 2023b.Citation: https://doi.org/10.5194/egusphere-2025-1828-RC2 -
AC1: 'Quick Reply on RC2', Bin Cao, 08 Sep 2025
reply
We thank Anonymous Referee #2 for the constructive comments. While waiting for other comments on our preprint before replying in detail, we want to briefly comment now on four issues:
(1) Anonymous Referee #2: "The application of FPM with lateral heat is provided in Sun et al. (2023)." This is very confusing to me for two reasons. First, vertical 1D models and cross sectional 2D models are usually very different types of models that implies different formulation of their physical equations, different numerical schemes for their resolution and different types of upper boundary conditions. The present study mentions surface energy balance calculation whereas Sun et al. (2023) forced their model with ground surface temperatures. So I do not understand how a 1D model can become a full 2D model (not a coupling of 1D simulation together).”
“Second, the given explanation is confusing. Sun et al. (2023) says "A 2D heat conduction model developed by Ling and Zhang (2004a) was used to simulate the permafrost thermal regime." So if it is the same model in 1D and 2D, is the present study actually presenting a new model or is it presenting improvements brought to a model published initially in 2004? Since the goal of this study is to present a new model, this kind of question need to clarified to consider publication.”
Response: We fully agree that 2D models are fundamentally different from 1D models. However, it is possible for a single model to support implementations in 1D, 2D, and even 3D. An example from the permafrost modeling community is the Control Volume Permafrost Model (CVPM) by Clow (2018), which “implements the nonlinear heat-transfer equations in 1-D, 2-D, and 3-D Cartesian coordinates, as well as in 1-D radial and 2-D cylindrical coordinates.”
We also acknowledge that referring to 2D capabilities may be misleading (Referee #1 raised a similar concern), as demonstrating model suitability requires additional applications. In the revised manuscript, we will remove the description of lateral heat transfer in FPM from the model description section and instead include it as part of the outlook.
(2) L62-104
“A physically-based surface energy balance scheme for different land surface cover types with varying snow regimes and properties was coupled to FPM, and was formulated as:
(1−α)Qsi +Qli +Qle +Qh +Qe +Qc = Qm”Major problem here. The study intends to describe the ground surface SEB (because we are in part 2.1 of a study presenting a permafrost model, with no mention to the snow scheme yet). Yet it looks like a description of a snow SEB, as evidenced by the expression “varying snow regime” and by the fact that the equation calculates energy for melt, which is inappropriate for a ground surface scheme. For a ground surface scheme, we want the SEB to give us access to the dE/dt of the surface so that it can force the heat diffusion/advection in the soil column. If we do not have that, we do not have the coupling between the climate and the temperature in the ground. This problem persist over the whole section 2.1. These are very important aspect of the model description that need to be carefully addressed so that we can understand what we are talking about.
Response: We agree that the term Qm is not appropriate here, as snow mass balance is not implemented in the current version of FPM. In the code, Qm is set to zero and therefore does not affect simulation results. In the revision, Qm will be removed.
Although both temperature and heat flux can serve as the upper boundary condition for soil heat conduction, FPM uses the surface temperature (Ts0)—whether snow or soil—similar to the approach in the CryoGrid community model (e.g., Eqs. 4 and 5 in Westermann et al., 2023). Accordingly, the surface energy balance will be revised as:
(1−α)Qsi +Qli +Qle +Qh +Qe = Qc (1)
and the sign of Qc should be revised,
Qc = (Ts0 – Tg)(zsn / ksn + zg / kg)-1 (10)
In FPM, Qle (Eq. 3) and Qh (Eq. 4) are functions of Ts0. Qe is derived as a function of Qn and Qc and thus also depends on Ts0. For each time step, Ts0—the upper boundary condition for the subsurface—is solved iteratively using the Newton–Raphson method to ensure energy conservation.
(3) L78 and 86
I have 2 problems with this quantification of the turbulent fluxes. First, if you take Preistley and Taylor for Q_e, then you have to take Q_h as the residual of the available energy for turbulent fluxes, which is:
Q_h=(1-α ∆/(∆+γ))×(Q_n-Q_g)
Otherwise your SEB is not energy conservative. Preistley and Taylor considers that the energy available for both turbulent fluxes is Q_n - Q_g and provide a formula to find what fraction of that available energy goes to the latent flux. So, in order to have a consistent surface energy balance budget, Q_h has to be the complementary fraction of Q_n - Q_g, not another formula based on another theory.
Response: We agree that using a consistent scheme for Qh and Qe would enhance model robustness. That said, combined approaches have been employed in previous studies. For instance, Song et al. used both Priestley–Taylor (for Qe) and Monin–Obukhov similarity theory (for Qh) in their surface energy balance model; other examples include Agam et al. (2010) and Kustas et al. (2003). Our detailed evaluations at both site and regional scales have demonstrated the suitability of FPM. We acknowledge that using two different theories may introduce additional uncertainties, and we will address this issue in the Discussion section in the revised manuscript.
Second, if you use Preistley and Taylor you will assume that the sum of the terms of your SEB equals 0 and not the energy variation of the soil surface. Otherwise you’d have to work with a modified version that would look like:
Q_e=α ∆/(∆+γ) (Q_n-Q_g-(∂E_surf)/∂t)
But then you would not be able to calculate together the turbulent flux and the energy variation of the surface (one too many unknown). Yet, it is through the energy variation of the surface that you can couple the SEB and the energy budget at depth in the ground (because the energy budget of the surface will drive the one of the subsurface). Therefore, this whole SEB description is very confusing to me and absolutely need to be fixed. For now I cannot understand the coupling with the subsurface.
Response: We agree that the current presentation of the SEB is misleading. As clarified above, energy conservation in the SEB is achieved through Eq. (1) and numerical solution via the Newton–Raphson iterative method.
(4) L43: I am surprised the authors do not mention other models used to study Tibetan permafrost like the GBEHM model that is greatly used within the Chinese community (Fang et al., 2025; Gao et al., 2018; Qin et al., 2017; Shi et al., 2020; Wang et al., 2023, 2018; Wang and Gao, 2022, 2025; Yang et al., 2023b, a), or the recent DHTC model (Linmao et al., 2024) and FLEXTopo-FS model (Gao et al., 2022).
Response: We agree that hydrological models applied in permafrost regions were not thoroughly reviewed in this section. This is primarily because such models typically include more detailed hydrological processes (e.g., water mass balance) but simplify the surface energy balance and heat conduction processes. For example, the DHTC model (Linmao et al., 2024) parameterizes ground heat conduction as a linear function of net radiation, and the FLEXTopo-FS model (Gao et al., 2022) uses the Stefan equation rather than a numerical solution for heat conduction.
Actually, Qin et al. (2017) used GIPL, not GBEHM. Nonetheless, both GBEHM (cited in L43 as Zheng et al., 2020) and GIPL (cited in L37 as Qin et al., 2017) are referenced in our manuscript.
References:
[1] Agam, N., Kustas, W. P., Anderson, M. C., Norman, J. M., Colaizzi, P. D., Howell, T. A., Prueger, J. H., Meyers, T. P., and Wilson, T. B. (2010): Application of the Priestley–Taylor Approach in a Two-Source Surface Energy Balance Model, Journal of Hydrometeorology, 11, 185–198, https://doi.org/10.1175/2009jhm1124.1.
[2] Clow, G. D. (2018): CVPM 1.1: a flexible heat-transfer modeling system for permafrost, Geosci. Model Dev., 11, 4889–4908, https://doi.org/10.5194/gmd-11-4889-2018.
[3] Gao, H., Han, C., Chen, R., Feng, Z., Wang, K., Fenicia, F., and Savenije, H. (2022): Frozen soil hydrological modeling for a mountainous catchment northeast of the Qinghai–Tibet Plateau, Hydrology and Earth System Sciences, 26, 4187–4208.
[4] Ling, F., and T. Zhang (2004), Modeling study of talik freeze-up and permafrost response under drained thaw lakes on the Alaskan Arctic Coastal Plain, J. Geophys. Res., 109, D01111.
Myhra, K. S., Westermann, S., and Etzelmüller, B. (2017) Modelled Distribution and Temporal Evolution of Permafrost in Steep Rock Walls Along a Latitudinal Transect in Norway by CryoGrid 2D. Permafrost and Periglac. Process., 28: 172–182.
[5] Kustas, W. P., Bindlish, R., French, A. N., and Schmugge, T. J.: Comparison of energy balance modeling schemes using microwave‐derived soil moisture and radiometric surface temperature, Water Resources Research, 39, 2003.
[6] Song, L., Kustas, W. P., Liu, S., Colaizzi, P. D., Nieto, H., Xu, Z., Ma, Y., Li, M., Xu, T., Agam, N., Tolk, J. A., and Evett, S. R. (2016): Applications of a thermal-based two-source energy balance model using Priestley-Taylor approach for surface temperature partitioning under advective conditions, Journal of Hydrology, 540, 574–587.Citation: https://doi.org/10.5194/egusphere-2025-1828-AC1
-
AC1: 'Quick Reply on RC2', Bin Cao, 08 Sep 2025
reply
Viewed
| HTML | XML | Total | Supplement | BibTeX | EndNote | |
|---|---|---|---|---|---|---|
| 543 | 62 | 17 | 622 | 18 | 10 | 27 |
- HTML: 543
- PDF: 62
- XML: 17
- Total: 622
- Supplement: 18
- BibTeX: 10
- EndNote: 27
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1