Bare Patches Created by Plateau Pikas Contribute to Warming Permafrost on the Tibet Plateau

Plateau pikas, small mammals native to the Qinghai‐Tibet Plateau (QTP), create bare patches through burrowing. No previous assessment exists on their impact on permafrost. This study fills this gap by simulating hypothetical scenarios in the Three Rivers Headwaters Region of the QTP using the Noah‐MP model for the plant growing seasons during 2015–2018. Our findings reveal a significant increase in soil temperature in the active layer due to pika‐induced bare patches, particularly during July–August. The average temperature rise at 2.5 cm depth was 0.36°C in permafrost regions and 0.29°C in seasonally frozen ground regions during August. Minimal impact on unfrozen water content was observed, with a slight increase in deep soil layers in permafrost regions, and negligible in seasonally frozen areas. These findings underscore the previously unexplored influence of pika burrowing on permafrost temperature, suggesting a potential risk of accelerating permafrost degradation, especially in permafrost‐dominated regions.


Introduction
Permafrost, defined as ground continuously frozen for two or more years (L.Zhao et al., 2010), forms a significant portion of the Qinghai-Tibet Plateau (QTP), accounting for approximately 67% of its total area (Cao et al., 2023).As the dominant vegetation type in the QTP, alpine grassland plays crucial roles in providing ecological functions such as biodiversity conservations, carbon storage, and permafrost preservation (Dong et al., 2020;Wen et al., 2013).Unfortunately, climate warming and human activities have contributed to the degradation of approximately 90% of alpine grassland, with 26% classified as extremely degraded (Harris, 2010;Shang & Long, 2007).
While climate change and overgrazing have significantly contributed to grassland degradation, recent studies highlight the role of rodent activities, particularly those of the plateau pika (Ochotona curzoniae), in grassland degradation (Qian et al., 2021;Zhang et al., 2020).As a keystone species in the Tibetan grassland, plateau pikas play a crucial role in enhancing plant species diversity and stabilizing vegetational communities (Smith & Foggin, 1999).Notably, the population of plateau pikas has increased in recent decades (Dong et al., 2013;Harris, 2010).Their subterranean lifestyle involves foraging and digging, bringing deep soil matters to the surface and creating distinct bare patches (Lai & Smith, 2003;Smith & Foggin, 1999).These bare patches alter radiation partitioning and ground heat conduction, affecting both thermal and hydrological conditions of the underlying permafrost (Qin et al., 2019;You et al., 2017).Previous plot-based studies indicate that pika-induced bare patches, with a positive correlation with pika burrow density, can elevate soil temperatures (Ma et al., 2018; J. Zhao et al., 2021).These findings suggest a potential contribution to permafrost degradation, although our understanding of this complex process remains limited.
Investigating the impact of pika-induced bare soil on permafrost hydrological-thermal processes at a regional scale requires accurate spatial distribution of these patches.However, their small size and spectral similarity to natural bare lands on satellite imagery pose challenges.To address this, we utilize a stochastic model to simulate pika burrow density, which informs the estimation of pika-induced bare patches based on statistical relationships derived from unmanned aerial vehicles (UAV) data.The Noah-MP land surface model (LSM) is then used to investigate the impacts of these bare patches on permafrost by configuring various scenarios in the Three Rivers Headwater Region (TRHR) of the QTP.This region is renowned for its unique role in China as an ecological shelter, extensive permafrost coverage, and widespread pika distribution (Yao et al., 2022).

Study Area and Data
The TRHR sits within the QTP hinterland, encompassing approximately 3.95 × 10 5 km 2 and bounded by longitude 89°24′-102°27′E and latitude 31°39′-36°10′E (Figure 1).Its elevation gradually declines from northwest to southeast, with an average of 4,000 m above sea level (a.s.l.) (Figure 1a).The mean annual air temperatures in the TRHR range from 5.6°C to 3.8°C, with a south-to-northwest decreasing trend (Yi et al., 2012).Spatially, annual precipitation varies significantly, with the southeast receiving 772 mm compared to 262 mm in the northwest.The majority of precipitation falls between June and September, influenced by warm and humid air currents emanating from the southern Bay of Bengal (Shi et al., 2017).The TRHR primarily consists of alpine steppe and meadow, although other cover types like bare land, wetland, shrub, and forest are also present (Figure 1b) (Chen et al., 2014).Permafrost distribution shows widespread coverage in the upper reaches of the Yangtze-River, whereas the eastern portion exhibits sporadic permafrost alongside seasonally frozen ground at lower elevations (Cao et al., 2023) (Figure 1c).The average active layer thickness in this region measured approximately 2.57 m between 1980 and 2010 (Xu & Wu, 2021).The growing season is relatively short, typically starting in early to late May and concluding in mid to late October (Dai et al., 2022).
The gridded Chinese Meteorological Forcing Data set (ITP-forcing) (0.1°in space and 3 hr in time) were used to drive the Noah-MP model (He et al., 2020).This data set includes seven meteorological elements including 2 m air temperature, 10 m wind speed, air pressure, specific humidity, precipitation rate, downward shortwave radiation, and longwave radiation.Vegetation type distribution was extracted from the 1:1 million scale Vegetation Atlas of China (Hou, 2019).Monthly leaf area index (LAI) data came from the level-4 Moderate Resolution Imaging Spectroradiometer (MODIS) 500 m global LAI and Fraction of Photosynthetically Active Radiation product (MCD15A2H).Stem area index (SAI) values, which are not directly obtainable from remote sensing, were sourced from the Noah-MP lookup table, considering 27 distinct vegetation types.Vegetation fractions were calculated from the normalized difference vegetation indices from the 1 km Terra MODIS vegetation indices data set (MOD13A3) (Ma et al., 2023), and bare land fractions were subsequently derived by subtracting vegetation fractions from 100%.A QTP soil data set with a 1-km resolution and 18 layers, encompassing the entire TRHR, was used (Wu & Nan, 2016).

Estimation of Pika-Induced Bare Patch Fractions
To estimate the fraction of pika-induced bare patches within each modeling cell, we first utilized a stochastic model to simulate pika burrow density across the study region.This model leverages information from a Bayesian additive regression trees (BART)-based species distribution model (Barbet-Massin et al., 2012;Carlson, 2020) that predicts pika habitat suitability.The habitat suitability map is categorized into four levels, unsuitable (0-0.2),low suitability (0.2-0.4), moderate suitability (0.4-0.6) and high suitability (0.6-1), based on predefined thresholds (Convertino et al., 2014;Mehdi & Arash, 2018).Three Gaussian distributions representing these suitability levels (low, moderate, high) were fitted based on observed burrow density data (Sun et al., 2015;Wei & Zhang, 2018).Unsuitable areas were considered to have no pika surface entrances.This stochastic model then generates pika burrow density for each grid cell.A more detailed introduction to this stochastic model can be found in Text S1 of Supporting Information S1.
Next, statistical relationships between pika burrow density and bare patch fraction were derived from UAV data (Zhang et al., 2021).Three linear equations were fitted for LBP, MBP and HBP levels (Figure S1 in Supporting Information S1).Each grid cell was categorized into LBP, MBP, or HBP based on satellite-derived bare land fractions from 2015 to 2018.Areas with NBP were merged into the LBP category in wake of satellite data's accuracy.By applying the corresponding UAV-based relationships to each grid cell based on its category, we estimated the pika-induced bare patch fraction.

The Noah-MP LSM
The Noah-MP model is distinguished by its incorporation of the latest parameterization schemes (Niu et al., 2011).It uses the Clapp-Hornberger water retention equation to solve vertical soil moisture distribution (Clapp & Hornberger, 1978).To represent land surface heterogeneity, it adopts a "semitile" subgrid scheme with separate energy budgets for vegetated and bare ground areas (F veg and 1-F veg , respectively).Equations S1-S3 in Supporting Information S1 detail these energy budgets.The model calculates net longwave radiation (L a ), latent heat (LE), sensible heat (H), and ground heat (G) fluxes (Equation S4 in Supporting Information S1) and enforces overall energy balance (S av + S ag = L a + LE + H + G).Vegetation canopy temperature (T v ), ground surface temperature in the vegetated fraction (T g,v ), and bare ground surface (T g,b ) are solved iteratively (Equations S1-S3 in Supporting Information S1).Noah-MP LSM allows users to choose among various parameterization options for individual processes, enabling customization based on research needs.The specific physical options used in this study are presented in Table S1 of Supporting Information S1.

Experiment Design
The presence of pika-induced bare patches affects three land surface parameters in Noah-MP: vegetation coverage, LAI, and SAI.Two scenarios were designed to assess the impact of bare patches on soil temperature and moisture dynamics.The first scenario (undisturbed vegetation) represents pre-disturbance conditions characterized by maximum vegetation coverage, LAI, and SAI.The process to determine maximum vegetation coverage includes selecting various environmental factors affecting vegetation growth, clustering grid cells based on these factors, and using the maximum vegetation fraction value within each cluster as a representative of the undisturbed state.LAI and SAI were then determined proportionally in relation to vegetation coverage, following the approach applied in Noah-MP.We considered 31 environmental factors, including terrain, temperature, and precipitation-related variables (details in Table S2 of Supporting Information S1) and applied principal component analysis (PCA) to these factors before clustering, utilizing the resulting PCA components for the clustering process.The Gaussian Mixture Model (GMM) clustering algorithm, known for its flexibility and ability to model various data distributions, was selected for spatial clustering (Reynolds, 2009).Cluster number was determined based on the automatic cluster detection feature of the GMM algorithm.
The second scenario (pika-induced bare patches) accounts for the impact of pika activity on vegetation cover.These pika-induced bare patch fractions were modeled based on statistical relationships and then subtracted from the pre-disturbance vegetation coverage to derive the post-disturbance state.Adjustments were made to both LAI and SAI in proportion to the vegetation fraction.To ensure realistic bare land proportions, the modeled pikainduced bare patches were constrained by satellite-derived bare land fractions.However, due to data limitations (Figure S1 in Supporting Information S1), the modeled pika-induced fractional bare patches only represent the average state for July-August from 2015 to 2018.To extend these values to the entire growing season (June through September), we linearly extrapolated the patch fractions for June and September using the ratio of the maximum vegetation fraction in that month relative to the mean vegetation fraction of July and August.The impacts were only considered during the plant growing seasons, as indicated by satellite-derived vegetation fraction analysis in the TRHR, which showed sharply reduced vegetation coverage outside these months.Consequently, during non-growing seasons, the exposed surfaces created by plateau pikas resemble unvegetated areas elsewhere, likely minimizing their influence on the underlying permafrost.

Model Settings, Calibration and Validation
The Noah-MP model was forced with the same ITP-forcing (0.1°) for both scenarios over the period of January 1984-December 2018.The model outputs during 2015-2018 were chosen for analysis, aligning with the timeframe of the UAV and in situ data.Prior to the start of the model run, a 300-year spin-up using the repeat forcing from 1979 to 1984 was conducted to mitigate the effects of initial values (Ji et al., 2022).The model was calibrated and validated using in situ soil temperature and moisture records from four permafrost monitoring sites (Figure 1b).Data from all years except the last served for calibration, with the final year reserved for validation.Six sensitive parameters (Cuntz et al., 2016;Hu et al., 2023): the exponent in Brooks-Corey relation (B parameter), soil porosity, saturated hydraulic conductivity, momentum roughness length, thermal conductivity of very fluffy snow, and albedo of fresh snow, were calibrated.All model parameters except for vegetation coverage, LAI and SAI, were held constant in both the undisturbed and pika-induced bare patch scenarios.The model's performance was evaluated using widely-used Nash-Sutcliffe efficiency coefficient (NSE) and Pearson's correlation (R).

Assessment of Model Performance
Noah-MP simulations of soil temperature and unfrozen water content were validated with independent in situ data from four permafrost sites (Figures S2 and S3 in Supporting Information S1).Overall, high R and NSE values for most depths and sites were observed, indicating the calibrated Noah-MP LSM successfully reproduced observed soil temperature and unfrozen water content.During the validation periods, soil temperature simulations exhibited high agreement with observations across all depths and sites, with R values exceeding 0.9 for all depth profiles.NSE values generally surpassed 0.68, except for deeper soil layers (>200 m) at certain sites (Figure S2l in Supporting Information S1).
The Noah-MP LSM effectively captured the temporal variations of unfrozen water content at all four sites and across different depths during the validation period.At the QT01 site, simulations demonstrated excellent agreement with observations across all depths, with R values ranging from 0.84 to 0.94 and NSE values ranging from 0.64 to 0.80 (Figures S3a-S3d in Supporting Information S1).While other sites displayed some discrepancies, the model still effectively captured the temporal trends.An outlier was observed at the QT08 site, 40 cm depth (Figure S3n in Supporting Information S1), where measured unfrozen water content remained consistently at 0. This potentially indicates a malfunction in the instrument, and therefore, R and NSE values were not calculated for this depth at the QTP08 site.

Spatial Distributions of Pika Burrows and Induced Bare Patches
Pika burrow density simulations revealed a distinct spatial pattern, with higher concentrations observed in the Yellow and Lancang River headwaters compared to the Yangtze River headwater (Figure 2a).Interestingly, the spatial distribution of pika-induced bare patches exhibited an inverse trend, with higher fractional bare patches concentrated in the Yangtze River headwater (Figures 2c and 2d).This contrasting pattern likely results from differences in climate conditions influencing vegetation recovery dynamics.The milder climate and abundant precipitation in the Yellow and Lancang River headwaters, while conducive to pika presence, also facilitate faster recovery of vegetation disturbed by burrow construction, leading to lower observed bare patch fractions.Conversely, the harsher climate and pre-existing bare lands in the Yangtze River headwater create conditions where minor pika disturbances can readily expand bare land areas, impeding vegetation recovery and resulting in higher observed bare patch fractions.Figure 2c depicts the 1-km resolution distribution of pikainduced bare patches, subsequently aggregated to 10 km to match Noah-MP simulations (Figure 2d).
Clustering analysis based on 31 environmental factors related to vegetation growth yielded five distinct clusters within the study region (Figure 2b).The attached table in Figure 2 presents the maximum vegetation coverage within each clsuters, which were late used to construct the undisturbed scenario.Clusters R1 and R2 displayed the highest vegetation cover, exceeding 0.94 in both July-August.Cluster R3 exhibited moderate cover, reaching 0.71-0.73during peak vegetation months.Moving westward, clusters R4 and R5 showed lower vegetation cover during peak months, ranging 0.27-0.50,with cluster R4 exhibiting slightly better conditions than cluster R5.

Impact of Pika-Induced Bare Patches on Soil Temperature and Moisture
Pika-induced bare patches significantly warmed shallow soil layers (2.5 and 66 cm), as shown in Figures 3a-3h and.Warming intensity varied spatially and temporally, with the most pronounced increases occurring in July-August.September also showed a stronger impact compared to June.Spatially, among the five vegetation clusters (Figure 2), cluster R4, with the most extensive pika-made bare patches, experienced the greatest warming, followed by R3.Conversely, cluster R5, naturally dominated by bare land, exhibited minimal change.In clusters R1 and R2 (except R2 north), pika bare patches had a weaker influence.At the regional scale, significant warming (>1°C) was concentrated in the south and southeastern Yangtze River headwater and the western Yellow River headwater, while minimal impact (<0.2°C) occurred in the western Yangtze River headwater and the southern Yellow River headwater.
These spatial patterns likely reflect the underlying permafrost and seasonally frozen ground distribution.Within permafrost zones, significant temperature increases were observed, with a maximum difference of 2.37°C between the contrasting scenarios and an average warming of 0.36°C at 2.5 cm depth in August (Figure 3c).In areas with seasonally frozen ground, notable differences were also observed in the northern Yellow River headwater, where the maximum difference reached 2.31°C and the average warming was 0.29°C, again at the shallow depth of 2.5 cm in August (Figure 3c).Notably, deeper layers (184 cm) showed minimal temperature changes, indicating limited impacts on the base of the active layer (Figures 3i-3l).
Plateau pika-induced bare patches had minimal impact on shallow soil moisture content (2.5 and 66 cm) (Figures 4a-4h).A slight increase in moisture content was observed at deeper layers (184 cm), nearing the bottom of the active layer (Figures 4j-4l).In space, the subtle changes were primarily concentrated in the southern Yangtze River headwater and the northern Yellow River headwater.The maximum observed increase was 0.06 m 3 m 3 (Figure 4k), with an average increase 0.007 m 3 m 3 in permafrost regions and 0.003 m 3 m 3 in seasonally frozen ground regions at 184 cm in August.These changes are considered negligible compared to overall soil moisture content.

Conclusions and Discussion
This study provides a new perspective for understanding permafrost responses to the bare patches caused by plateau pikas over the QTP through hypothetical numerical experiments.The following conclusions have been drawn.
1.While plateau pikas prefer warmer and wetter headwater regions of the Yellow and Lancang Rivers, pikainduced bare patches are more extensive in the drier and colder Yangtze River headwater.This inversely related pattern reflects the complex interplay between various environmental conditions and vegetation recovery across these regions.2. The pika-created bare patches led to a substantial soil temperature rise, especially in shallow soil layers (2.5 and 66 cm).The impact was greater in permafrost zones compared to seasonally frozen ground.In August, the average temperature rise at 2.5 cm in the permafrost zone was 0.36°C (2015-2018), compared to 0.29°C in the seasonally frozen ground zone.

Geophysical Research Letters
10.1029/2024GL108976 3. The influence of pika-induced bare patches on soil unfrozen water was relatively minor.A slight increase in moisture content was observed near the active layer base, averaging 0.007 m 3 m 3 at 184 cm in August simulations.
Uncertainties exists in estimating pika burrow density using the stochastic model.More field surveys could improve these estimates.Uncertainties are also associated with estimating pika-induced bare patch fractions and determining undistributed vegetation conditions.Rising permafrost temperature due to pika bare patches, combined with climate warming, would further amplify the risk of permafrost degradation.Pika burrow systems, averaging 55.2 cm deep (Wei & Zhang, 2018), act as macropores in the active layer, facilitating deeper water infiltration.These macropores potentially create pathways for convective heat transfer, allowing warm precipitation to reach deep soil layers and accelerate permafrost degradation (Luo et al., 2020;Qin et al., 2019;Wilson & Smith, 2015).

Figure 1 .
Figure 1.Maps showing the location, topography (a), land cover types (b), and frozen ground types (c) in the Three River Headwaters Region on the Qinghai-Tibet Plateau.Labels 1, 2 and 3 in (a) indicate the headwaters of the Yangtze River, Lancang River, and Yellow River, respectively.Pika presence locations were recorded based on previous studies spanning 2011-2018.The land cover map is adapted from Chen et al. (2014), and the frozen ground map from Cao et al. (2023).

Figure 2 .
Figure 2. Modeled pika burrow density distribution and induced bare patch fractions in the Three Rivers Headwaters Region (TRHR).(a) Spatial distribution of pika burrows across the TRHR, (b) five clusters based on environmental factors used to derive undisturbed vegetation fractions, (c) modeled pika-induced fractional bare patches at 1-km resolution, and (d) Pika-induced fractional bare patches aggregated to 10 km from the 1-km map (c)."No data" areas in (c) and (d) indicate unsuitable regions for pika habitat and were excluded from simulations of pika-induced bare patches.The accompanying table presents the maximum vegetation fractions for each cluster (b) during the growing season (June-September).These values were used to construct the undisturbed scenario.

Figure 3 .
Figure 3. Maps illustrating changes in soil temperature during the plant growing seasons induced by the presence of pika-induced bare patches.Each row represents a specific soil depth.Values were computed by subtracting soil temperatures simulated under pika-induced bare patch scenario from those simulated under undisturbed vegetation scenarios.Positive values indicate areas where pika-induced bare patches led to increase in soil temperature, while negative values represent areas where they caused decreases.

Figure 4 .
Figure 4. Maps illustrating differences in soil liquid water content (SWC) during the plant growing seasons induced by the presence of pika-induced bare patches.Each row represents a specific soil depth.Values were computed by subtracting SWC simulated under pika-induced bare patch scenario from those simulated under undisturbed vegetation scenarios.Positive values indicate areas where pika-induced bare patches led to increase in SWC, while negative values represent areas where they caused decreases.