Quantifying permafrost ground ice contents in the Tien Shan and Pamir (Central Asia): a petrophysical joint inversion approach using a geometric mean model

In the Central Asian Tien Shan and Pamir mountain ranges, permafrost is extensive, but in situ data on permafrost remain scarce. Quantitative analysis of permafrost's subsurface components – ice, water, air, and rock – is vital for not only discerning the impact of climate change on increased slope instability due to permafrost degradation, but also understanding its role as a potential water resource in high-altitude environments. Recent studies have employed a petrophysical joint inversion (PJI) approach combining geoelectrical and seismic refraction data to model the subsurface's four phases (fractions of air, water, ice, and rock). However, most of these studies primarily rely on Archie’s law, which has limitations in coarse blocky substrates typical of mountainous terrains. Recognizing this limitation, the electrical geometric mean (PJI-GM) model may be used as an alternative implementation within PJI. In this study, we assess the suitability of using the PJI-GM model across an extensive geophysical dataset comprising 22 profiles in Central Asia (Kyrgyzstan and Tajikistan). Our goals are to (i) address the existing data gap concerning mountain permafrost and ground ice contents in the Tien Shan and Pamir of Central Asia and (ii) evaluate the performance of the PJI-GM model in comparison to Archie's law within the PJI framework across the different landforms at remote sites. The findings reveal that the ground ice content is more specific to landform types than to the different geographic regions surveyed, with rock glaciers exhibiting the highest mean ice contents (38 %–60 %), followed by moraines (18 %–40 %), talus slopes (20 %–40 %), and fine-grained sediments (0 %–20 %). The PJI-GM model performed especially well for ice-rich landforms such as rock glaciers, accurately reflecting high ice contents with minimal variability between different model runs. The quality of a model result was here assessed by comparing a multitude of different model runs with different sets of inversion parameters and petrophysical variables using a clustering approach. This research provides one of the first comprehensive (geophysical) in situ datasets on permafrost on various landforms and sites in Central Asia, highlighting the potential of the PJI-GM model as a more suitable alternative to Archie’s law, particularly for rock glaciers and other ice-rich landforms. These findings significantly advance our understanding of permafrost in the Tien Shan and Pamir and serve as a baseline dataset for future modeling studies.

1 Introduction

With ongoing climate change, permafrost – defined as subsurface material with temperatures at or below 0 °C for at least 2 consecutive years – is experiencing increased warming and degradation (Biskaborn et al., 2019; Etzelmüller et al., 2020). In mountain regions, this degradation has two key impacts. First, it affects slope stability, compromising mountain slopes and increasing risks of landslides, debris flows, and altered sediment transport (Daanen et al., 2011; Haeberli et al., 2017; Ravanel et al., 2017). Second, it raises questions about its role in the hydrological cycle, including its contribution as a water resource and its influence on evaporation and runoff patterns (e.g., Jones et al., 2019; Luo et al., 2020; Martin et al., 2023). Although permafrost may offer an alternative water resource, especially in dry regions, its role in seasonal water supply and river system sustainability remains unclear (Arenson and Jakob, 2010; Arenson et al., 2022; Amschwand et al., 2024). Comprehensive data on permafrost distribution, thermal state, and ground ice content are crucial for understanding climate impacts on slope stability and hydrology, as well as for mitigating geohazards. In the Tien Shan and Pamir mountain ranges of Central Asia, despite the expected widespread occurrence of permafrost (Marchenko et al., 2007), a considerable data gap exists (Barandun et al., 2020; Hoelzle et al., 2019). This scarcity of data underscores the urgent need for comprehensive in situ investigations to characterize permafrost distribution, thermal state, and ground ice content in the Tien Shan and Pamir mountain ranges.

In situ methods are crucial for effectively detecting and monitoring changes in the thermal state and ground ice contents of mountain permafrost, as most permafrost features are not directly detectable with remote sensing techniques. Permafrost is typically monitored using a combination of methods such as boreholes, geophysical measurements, and ground surface temperature (GST) loggers (e.g., PERMOS, 2023). Boreholes enable direct measurements of ground temperatures for thermal analysis, and drill cores can provide information about ground ice contents, thereby yielding accurate data on the subsurface characteristics of permafrost (Noetzli et al., 2021). In Central Asia, despite limited current data, early permafrost research was extensive, with systematic studies in the northern and central Tien Shan of Kyrgyzstan and Kazakhstan beginning in the 1950s (e.g., Ermolin et al., 1989; Gorbunov, 1967, 1970), but Duishonakunov (2014) noted that there were permafrost observations in this region even earlier. However, much of this research diminished after the 1990s, resulting in significant gaps in long-term monitoring efforts. Until the 1990s, numerous boreholes were being drilled in this region (Marchenko et al., 2007). A lower boundary for sporadic permafrost in the Tien Shan was estimated to be at an altitude of 2800–3000 m a.s.l. (Gorbunov et al., 1996). Marchenko et al. (2007) found an increase in permafrost temperatures in the Tien Shan within the range of 0.3 to 0.6 °C using observations with an average increase in the active layer by about 23 % from the 1970s to 2004. Similarly, Seversky (2017) detected a slight increasing trend in ground temperatures between 1995 and 2016 (0.01 °C a⁻¹ from 1974 to 2016) in boreholes in the northern Tien Shan. While this historical basis provides valuable context, accessing the datasets can be challenging and only one borehole in the northern Tien Shan is still actively observed today within the GTN-P network (Seversky, 2017). Furthermore, borehole drilling is expensive, is logistically challenging, and can only provide point information.

Geophysical techniques such as electrical resistivity tomography (ERT), refraction seismic tomography (RST), or ground-penetrating radar (GPR) are invaluable for expanding our understanding of permafrost over larger areas. They can, for example, be used to delineate the active layer and taliks, as well as assess changes in subsurface properties if measurements are repeated regularly (e.g., Hauck and Kneisel, 2008; Hilbich et al., 2009; Monnier and Kinnard, 2013; Mollaret et al., 2019; Halla et al., 2021; Mollaret et al., 2020; Kneisel et al., 2008; Vonder Mühll et al., 2001; Boaga et al., 2020; Herring et al., 2023). In Central Asia, geophysical investigations are limited to the central and northern Tien Shan (Seversky, 2017). However, spatial and temporal coverage is relatively sparse. ERT surveys were conducted in the northern Tien Shan in 2013 and in 2017, where high resistivities typical of (ice-rich) permafrost were found. Bolch et al. (2019) used GPR measurements to estimate ground ice contents of ice–debris complexes in the central Tien Shan and mapped a total of 74 rock glaciers and ice–debris complexes using remote sensing. Boreholes have revealed the presence of ground ice in the form of ice lenses in moraines (Marchenko et al., 2007). Bolch and Marchenko (2009) estimated rock glacier ground ice contents in the Tien Shan based on an empirical relationship proposed by Brenning (2005). However, quantitative estimates of ground ice contents based on field measurements are, to our knowledge, not yet available for the Tien Shan. Furthermore, most remote sensing research on permafrost in the area is focused on rock glaciers and is often lacking in situ data for validation (Blöthe et al., 2019; Sorg et al., 2015; Kääb et al., 2021; RGIK, 2023).

Data and research on permafrost in general are much scarcer in the Pamir and the Pamir Alay (Barandun et al., 2020). Here, permafrost occurrence has been described down to an altitude of 3800 m a.s.l. Climate model projections for the extended Tibetan Plateau, which in most studies includes parts of the Pamir, suggest a reduction in near-surface permafrost of 39 % by 2050 and up to 81 % by 2100 (Bolch et al., 2019). Most current studies that focus on permafrost in the Pamir rely on remote sensing and focus not on ground ice contents but rather on hazards associated with potentially increasing slope instability as permafrost degrades (Jones et al., 2021b; Mergili et al., 2013; Mergili and Schneider, 2011). To our best knowledge, no prior geophysical investigations on permafrost have been conducted in the Pamir region to date.

Yet, understanding permafrost subsurface conditions, such as volumetric contents of ice, water, and rock, is crucial for developing a comprehensive understanding of permafrost processes and for evaluating associated degradation risks. Permafrost genesis in mountainous regions is influenced by a combination of climatic, geological, and environmental factors (Gilbert et al., 2016). The formation of ground ice within permafrost can occur through processes such as the freezing of infiltrating precipitation, the migration of water to freezing fronts, the burial of glacier ice, and the burial of snow (e.g., Pollard, 1990; Bockheim and Tarnocai, 1998; Monnier et al., 2011; Kenner et al., 2017; Gilbert et al., 2016). These processes vary across different landforms, leading to diverse permafrost characteristics and ground ice contents. However, quantifying ground ice content within different landforms remains a significant challenge. For example, rock glaciers have been shown in numerous studies to contain significant amounts of ground ice (e.g., Vonder Mühll and Holub, 1992; Vonder Muhll and Haeberli, 1990; Krainer et al., 2015; Monnier et al., 2011). However, even for relatively well studied landforms like rock glaciers, most estimates rely on empirical relationships to estimate ground ice contents (Brenning, 2005), with only a limited number of studies providing quantitative measurements (e.g., Halla et al., 2021; Pavoni et al., 2023; Hilbich et al., 2022; Mollaret et al., 2020). This lack of direct measurement data and the complexity of permafrost subsurface conditions make it difficult to establish a consensus on typical ice content ranges, particularly for landforms beyond rock glaciers such as talus slopes or moraines. To our knowledge, no single reference paper comprehensively summarizes expected ground ice contents across various mountain permafrost landforms.

However, in recent years, various (geophysical) approaches have been employed to quantify these subsurface constituents, in particular the ground ice contents. Hauck et al. (2011) introduced the Four Phase Model (4PM), which integrates electrical resistivity and P-wave velocity measurements to characterize water, ice, air, and rock contents in the subsurface. This model has been used in various studies, and results have been shown to fit well with ground truth data such as borehole data or field observations (e.g., Halla et al., 2021; Hilbich et al., 2022; Kunz et al., 2022). The 4PM uses Archie’s law (Archie, 1942) to relate the bulk resistivity of a material to its porosity and the resistivity of the pore water, which works well in environments where electrolytic conduction dominates. However, Archie’s law does not directly account for the fractions of air and ice. In the 4PM, these fractions are indirectly constrained by integrating ERT data with RST data, which provides additional information on the seismic properties of the subsurface materials. Wagner et al. (2019) further developed a petrophysical joint inversion (PJI) scheme that builds on the principles of the 4PM but jointly inverts the electrical and seismic data. The PJI scheme has been applied in multiple studies (e.g., Mollaret et al., 2020; de Pasquale et al., 2022; Pavoni et al., 2023; Steiner et al., 2021b; Klahold et al., 2021). Mollaret et al. (2020) tested different petrophysical models within the PJI scheme and suggested that the so-called electric geometric mean model (hereafter referred to as PJI-GM) could offer more realistic results compared to the commonly used Archie's law, which is currently the main resistivity equation implemented in the PJI scheme. This stems from the recognition that Archie's law (hereafter referred to as PJI-AR) is generally considered valid when electrical conduction through fluids within the pore space dominates over conduction through the solid matrix itself, a condition that is not universally justified. The geometric mean model assumes that the modeled space is composed of a mixture of the four phases (rock, ice, water, and air), with each phase being randomly distributed within the subsurface. This approach allows for the inclusion of the fractions of all four phases, as opposed to Archie's law. However, the PJI-GM model presents challenges as additional unknowns (the resistivities of ice, air, and rock in addition to the pore water resistivity that is needed in Archie's law) are added to the system of equations. Furthermore, securing model convergence within PJI-GM can be challenging, as indicated by substantial errors in numerous model outputs (Mollaret et al., 2020). Finally, PJI-GM has not yet been tested extensively for a large dataset comprising multiple examples of different landforms.

Since 2021, we have been addressing the in situ data gap on permafrost in Central Asia by conducting extensive geophysical surveys. These surveys involved ERT and RST measurements at various study sites in the Tien Shan and Pamir mountain ranges. In this study, we present these data and assess the suitability of using the PJI-GM model within the PJI framework to estimate ground ice content distribution across an extensive geophysical dataset comprising 22 profiles in Central Asia (Kyrgyzstan and Tajikistan). Our research encompasses diverse landforms, including moraines, rock glaciers, talus slopes, and fine-grained sediments. Our goals are to (i) address the existing data gap concerning mountain permafrost and ground ice contents in the Central Asian region and (ii) evaluate the performance of the geometric mean model in comparison to Archie's law across different landforms. The baseline dataset established within this study is essential for developing accurate models and predictions of future permafrost dynamics and their associated impacts on hydrology and geohazards in the face of climate change and may assist local populations in adapting to forthcoming changes by informing water resource management and disaster risk reduction policies.

2 Study sites

Geophysical measurements were carried out at 10 different sites distributed across the Pamir and Tien Shan mountain ranges and ranging from 3100 to 4580 m a.s.l. A total of 38 ERT profiles and 22 RST profiles were measured on different landforms, which were categorized into rock glaciers (RGs), talus slopes (TSs), fine-grained sediments (SEDs), and moraines (MOs). Figure 1 shows the location of each study site with pictures of the different landforms investigated. The climate in the Central Asian region is mostly semi-arid, with high seasonal variability due to its continental location and considerable regional variations (e.g., Aizen et al., 2009; Haag et al., 2019; Barandun and Pohl, 2023). The sites chosen for the permafrost analyses reflect different climatic and geomorphological settings. They are part of a comprehensive cryospheric monitoring network being established in the region, covering all cryospheric variables (snow, glaciers, permafrost). Meteorological stations were installed near most of the study sites as part of projects to re-establish the glacier monitoring network from past and ongoing projects efforts (Hoelzle et al., 2017; Schöne et al., 2013; Zech et al., 2021). Mean annual air temperatures (MAATs) and mean annual precipitation were calculated. The results are summarized in Table 1, which provides a summary of the study sites, highlighting variations in climate and other relevant site information. A short introduction to each sub-region (after Barandun and Pohl (2023) and Zandler et al. (2016)) and corresponding study sites is given in the following sections. Figure 2 shows a zoomed-in view of each study site, showing the location of all geophysical profiles (ERT and RST). Data acquisition details for each geophysical profile are given in Table 2.

Figure 1 Study sites are marked with pink circles, with diameters proportional to the number of geophysical profiles. Each circle is labeled with its study site name and includes the number of ERT and RST profiles conducted. Different sub-regions are distinguished by different colors. The pictures display sampled landforms such as talus slopes, moraines, rock glaciers, and fine-grained sediments. Sources of the background maps: Esri, DigitalGlobe, GeoEye, i-cubed, USDA FSA, USGS, AEX, Getmapping, Aerogrid, IGN, IGP, swisstopo, and the GIS User Community. Publisher's remark: please note that the above figure contains disputed territories.

2.1 Northern Tien Shan

The northern and northwestern Tien Shan sub-region is characterized by a comparatively moist climate with precipitation rates of about 700 mm a⁻¹ (e.g., Aizen et al., 2006; Barandun et al., 2018; Guan et al., 2022). Earlier studies have identified a significant number of large rock glaciers in the region (Blöthe et al., 2019; Bertone et al., 2019). The study sites Golubin and no. 599 are located within this sub-region. MAAT at the Golubin site (2014–2021, at an altitude of 3305 m a.s.l.) is approximately −1.47 °C (see Fig. 1). The profiles located at the Golubin study site in the Ala Archa catchment (about 30 km from the capital of Kyrgyzstan, Bishkek) range from 3050 to 3410 m a.s.l. and include all four landform categories. Study site no. 599 is located to the north of Lake Issykul. The ERT profile measured at this site is located on a moraine at an altitude of 3780 m a.s.l. Unfortunately, no meteorological station is available in proximity to this site.

Hoelzle et al. (2017)Kronenberg et al. (2020)Schöne et al. (2013)Aizen et al. (2006)Machguth et al. (2024)Machguth et al. (2024)Hoelzle et al. (2017)Hoelzle et al. (2017)

Table 1 Study site overview: MAAT (mean annual air temperature) and mean annual precipitation values were derived from meteorological station (MS) data located near the geophysical measurement sites. These stations are described in Hoelzle et al. (2017) and Schöne et al. (2013). The time span for the calculation is indicated in parentheses next to the MS name. Wherever this was not possible, the information was taken from other studies and is indicated in the “References” column.

2.2 Central Tien Shan

Compared to the northern Tien Shan, the climate in the central Tien Shan is drier with annual precipitation amounts of about 350 mm a⁻¹ (Aizen et al., 2006). The study sites Kumtor and Suek are located within this sub-region south of Lake Issykul. Both study sites are located on a high-mountain plateau (mean of 3600 m a.s.l.) in the Upper Naryn catchment. The profiles measured at the Kumtor site are all located on fine-grained, vegetated sediments, about 5 km from the Kumtor gold mine. Several boreholes were drilled at this location in the 1980s (Marchenko et al., 2007) but are mostly inactive today. A new borehole was drilled in 2022 at the location of profile KUM04, revealing frozen conditions, a shallow active layer of approximately 1–1.5 m, and saturated ground ice conditions in the upper part of the drill core. The profiles measured at the Suek study site cover all four landform types and are distributed along the Suek pass with altitudes ranging from 3400–3880 m a.s.l.

Figure 2 Distribution of geophysical profile lines (spatial extent and orientation) across each study site. Each line is labeled with the corresponding profile name, which can be cross-referenced with the information provided in Table 2. The locations of the meteorological station (MS) (mentioned in Table 1) are marked by stars (pink and yellow). If the station is located more than 10 km away from the profiles, the direction and distance to the station are indicated with the arrows above the yellow stars. Sources of the background maps: Esri, DigitalGlobe, GeoEye, i-cubed, USDA FSA, USGS, AEX, Getmapping, Aerogrid, IGN, IGP, swisstopo, and the GIS User Community.

Table 2 Data acquisition parameters and on-site permafrost observations. The profiles which have both ERT and RST data are marked in bold. RST profiles are usually shorter than the ERT profiles, but we used the same spacing for both methods. The lines are shown in Fig. 2. ERT configurations: W, Wenner; DD, dipole–dipole. Landform classes: RG, rock glacier; MO, moraine; TS, talus slope; SED, fine-grained sediment. Only the Kumtor site has previous borehole data that were available for comparison for this study.

2.3 Pamir Alay

The Pamir Alay meridionally separates the Pamir and Tien Shan mountains. It is located in the northwest of the Pamir mountain range and encompasses different mountain ranges spanning southern Kyrgyzstan to the northwest of Tajikistan. The study sites Abramov (KG) and Yakarcha (TJ) are located in this sub-region. At the Abramov study site, located in the Vakhsh catchment, MAAT at an altitude of 4100 m a.s.l. is −4.6 °C (measurements from 2010–2021). Typical precipitation rates in this part of the Pamir Alay are around 750 mm a⁻¹ (Barandun et al., 2015; Kronenberg et al., 2020). The geophysical surveys include various rock glaciers, talus slopes, and fine-grained sediments, as well as multiple measurements on the Abramov glacier lateral moraine. At the Yakarcha study site, located to the north of Dushanbe in the western part of the Pamir Alay in the Varzob catchment, MAAT is around −2.14 °C (measurements from 2019–2023) at an altitude of around 3500 m a.s.l. Precipitation in the region (measured at the nearby Ansob pass) is around 400 mm a⁻¹ (Rahmonov et al., 2017). Here, ERT profiles were measured on a large rock glacier.

2.4 Western Pamir

The western Pamir is characterized by deeply incised valleys and high mountain ranges (5000–6000 m a.s.l.) (Breu et al., 2003). The climate is mainly influenced by the westerlies with minimal rainfall in the summer months (Aizen et al., 2009; Barandun and Pohl, 2023). The western Pamir shows extreme precipitation differences at regional to local scales, e.g., with the highest monitored long-term mean of 2234 mm a⁻¹ at Fedchenko Glacier (4300 m a.s.l.) in comparison to the Lake Sarez station 50 km to the south (3290 m a.s.l.), where average precipitation is only about 110 mm a⁻¹. Large glaciers (e.g., Fedchenko) and rock glaciers are abundant in the western Pamir. The study site no. 457 is located within this sub-region. Geophysical profiles are located on a rock glacier and on fine-grained, mostly vegetated sediment at a mean altitude of 4580 m a.s.l. The closest meteorological station is located in the village of Jelondy, about 30 km away from the site. Precipitation is not measured at this station. MAAT in the years 2020–2022 at an altitude of 3560 m a.s.l. was −1.11 °C.

2.5 Eastern Pamir

The westerlies in combination with the Indian summer monsoon (ISM) are the main drivers of the climate in the eastern Pamir (Aizen et al., 2009). There is a negative west–east precipitation gradient (Fuchs et al., 2013; Pohl et al., 2015), making the eastern Pamir the driest of all the regions investigated in this study, with very little overall precipitation (40–140 mm a⁻¹) (Breu et al., 2003; Pohl et al., 2015; Barandun and Pohl, 2023). Most of the eastern Pamir is further characterized by a high plateau (mean of 4000 m a.s.l.) with wide valleys. While large rock glaciers are present in the southeastern Pamir, towards the northeastern Pamir they diminish in size and frequency. The study sites Bulunkul, Zulmart, and Karakul are located in this sub-region, distributed along a north–south axis. Bulunkul in general is a special site due to geographic barriers surrounding the site, resulting in low precipitation and exceptionally low temperatures (Pohl et al., 2015). The geophysical profile (at an altitude of 3720 m a.s.l.) is located on fine-grained, partly vegetated sediment. The profiles at the Zulmart study site are located on rock glaciers, fine-grained sediment, and a moraine, all at a similar altitude of about 4530 m a.s.l., with mean temperatures of about −4.1 °C (measured since 2021). The Karakul study site is located to the north of Lake Karakul at an altitude of 4200 m a.s.l. Measurements (1934–2017) from the meteorological station in the village of Karakul, located 20 km to the south of the study site at an altitude of 3900 m a.s.l., indicate MAAT of −3.7 °C and low precipitation of about 80 mm a⁻¹. The geophysical profiles are located on fine-grained sediments.

3 Methods

We use two well-established geophysical techniques, ERT and RST, and a modified PJI approach to assess the permafrost distribution and the ground ice contents at the different study sites. The steps of our methodology are summarized in Fig. 3 and will be introduced in the following sections.

Figure 3 Method flowchart showing the different steps of the methodological approach presented in this paper to estimate air, water, ice, and rock contents of different permafrost landforms.

3.1 Electrical resistivity tomography (ERT) and refraction seismic tomography (RST)

The ERT data were collected using a Syscal instrument (Iris instruments) with electrode spacing ranging from 2 to 5 m, primarily employing the Wenner configuration for its good signal-to-noise ratio while also carrying out supplementary dipole–dipole measurements for enhanced lateral resolution at some of the profiles (see Table 2). To efficiently compare the gathered datasets from diverse locations and landforms, it is essential to ensure consistent processing and analysis of the data (Herring et al., 2023; Mollaret et al., 2019). Therefore, the measured apparent resistivity data were filtered with criteria proposed by Mollaret et al. (2019) to eliminate outliers and to avoid overfitting of random noise contained in the data. Despite partly challenging terrain conditions, more than 90 % of the original quadrupoles remained after this filtering approach, indicating good overall quality of the raw data. A summary of the filtering statistics is presented in Table A1. The filtering is followed by iterative data inversion using the open-source pyGIMLi library with a smoothness-constrained least-squares generalized Gauss–Newton algorithm (Rücker et al., 2017). Regularization parameters for individual inversions include the smoothness regularization parameter (α) and a parameter for the relative weight for vertical boundaries (zWeight). These parameters were individually selected for all profiles through a sensitivity analysis and using the L-curve method proposed in Rücker et al. (2017) to optimize model response. The evaluation of inversion quality involved assessing the dimensionless error-weighted χ² parameter, which quantifies the misfit between the model response and the data for a given data error, along with the RMSE, which provides a measure of average magnitude of errors between observed data and model response. The hereby optimized regularization parameters α and zWeight used in the individually inverted ERT data are in the following kept consistent with those used in the PJI runs.

Co-located on 22 of the 38 ERT profiles, we conducted refraction seismic tomography (RST) surveys using a geode system equipped with 24 geophones, also spaced 2 to 5 m apart. RST first-arrival picking was performed using the software ReflexW (Sandmeier, 2024). Picking the first arrivals for each geophone can be challenging when dealing with poor data quality, which may result from inadequate anchorage of the geophones in the ground or other disruptive factors such as strong winds causing noise. Nevertheless, the correct identification of first-arrival travel times is a critical step in RST data quality. To ensure good quality, only datasets with over 80 % confidently identified first arrivals were considered suitable for further processing and subsequent data inversion to ensure reliable results. Only one RST dataset had to be excluded from the further processing steps because of bad data quality. Similarly to ERT, data inversion was carried out within pyGIMLi, where the regularization parameters were chosen in the same way as for the ERT inversions, explained above. The quality of the inversion results was assessed by forward modeling of the ray paths and subsequent comparison with measured travel times. A summary of the RST filtering and inversion statistics can be found in Table B1. It has to be noted that RST surveys take longer to conduct compared to ERT surveys, so they were performed only at specific ERT profiles as shown in Fig. 2, and they typically do not cover the entire length of the ERT profiles.

3.2 Petrophysical joint inversion

To quantify ground ice content, we employ the petrophysical joint inversion (PJI) approach, as developed by Wagner et al. (2019). The PJI model combines a set of petrophysical equations to quantify subsurface ice, water, and air content based on measured seismic P-wave travel times and electrical resistivities. While porosity is a necessary input in the original 4PM formulation (Hauck et al., 2011), it is often poorly known. Furthermore, utilizing independently inverted seismic and electrical data can yield non-physical results (Wagner et al., 2019; Mollaret et al., 2020). The PJI model offers the advantage of a simultaneous and physically consistent inversion for rock, ice, water, and air contents using the apparent resistivities and travel times as input data. The underlying assumption in the model is that the subsurface is composed of four phases: rock matrix (f_r), water (f_w), air (f_a), and ice (f_i):

$$\begin{array}{}\text{(1)}& f_{\mathrm{w}}+f_{\mathrm{r}}+f_{\mathrm{i}}+f_{\mathrm{a}}=\mathrm{1}.\end{array}$$

We employ Eq. (2), known as the time-average equation, to estimate the bulk velocity v based on the constituent fractions and their respective velocities. This method is an expansion of Timur's time-averaging approach (Timur, 1968) to accommodate all four phases found in permafrost (Hauck et al., 2011). Equation (2) expresses the notion that the inverse of the P-wave velocity $\mathrm{1}/v$ (slowness) in a mixture is equivalent to the combined slownesses of each component, weighted by their respective volumetric fractions:

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\mathrm{1}}{v}}={\displaystyle \frac{f_{\mathrm{w}}}{v_{\mathrm{w}}}}+{\displaystyle \frac{f_{\mathrm{r}}}{v_{\mathrm{r}}}}+{\displaystyle \frac{f_{\mathrm{i}}}{v_{\mathrm{i}}}}+{\displaystyle \frac{f_{\mathrm{a}}}{v_{\mathrm{a}}}}.\end{array}$$

The phase velocities v_w (1500 m s⁻¹), v_i (3500 m s⁻¹), and v_a (300 m s⁻¹) were considered constant and equal to values that are well established in literature (e.g., Hauck and Kneisel, 2008). However, v_r is site-dependent and could only be defined with laboratory tests, which are not available at the sites investigated in this study. Therefore, we also considered v_r constant for all sites, with an average value of 5500 m s⁻¹.

3.2.1 Petrophysical equations for electrical resistivity

Archie's law

Archie's law (Archie, 1942) is the most commonly used petrophysical equation; it relates the bulk resistivity ρ to the pore water resistivity ρ_w, the porosity, and the fraction of the pore space occupied by liquid water:

$$\begin{array}{}\text{(3)}& \mathit{\rho }={\mathit{\rho }}_{\mathrm{w}}(\mathrm{1}-f_{\mathrm{r}}{)}^{-m}{\left({\displaystyle \frac{f_{\mathrm{w}}}{\mathrm{1}-f_{\mathrm{r}}}}\right)}^{-n},\end{array}$$

where ρ_w is resistivity of pore water, n is the saturation exponent, and m is the cementation exponent. The exponents m and n are substrate-specific parameters that are assumed to be constant over space in our study, as no detailed subsurface information is available, which is generally the case for the study sites in Central Asia. However, Archie's law is considered valid only when electrolytic conduction dominates. This is not universally justified for different substrates and landforms in mountainous terrain, especially for coarse blocky material (Duvillard et al., 2018; Coperey et al., 2019). Archie's law relies on the assumption that electrical conduction occurs through the water in the pore space, which can lead to an underestimation of porosity when non-conductive phases, such as ice or air, are dominant. Also, other conduction mechanisms such as surface conduction in an electrical double layer or at the interface between the ice and the unfrozen material are neglected (Maierhofer et al., 2024). The advantages of Archie’s law include its simplicity and long-established application in hydrogeology, but it can lead to misinterpretations of resistivity in mountain permafrost environments, often resulting in an overestimation of rock content, due to its inability to account for non-conductive materials. Furthermore, the fractions of ice and air are not included in Archie's law (see Eq. 3) and are therefore not constrained by the equation. The PJI version using Archie's law is hereafter called PJI-AR.

Geometric mean model

The PJI model using the electrical geometric mean (PJI-GM) offers an alternative petrophysical equation to link the measured bulk resistivity to the volumetric fractions of the four phases. It has the advantage of including rock, ice, and air resistivities in addition to the resistivity of the pore water resistivity (${\mathit{\rho }}_{\mathrm{r},\mathrm{i},\mathrm{w},\mathrm{a}}$ in Eq. 4); (Mollaret et al., 2020), similarly to the P-wave velocity equation. It represents an alternative approach that might be better suited for environments with mixed conductive and non-conductive phases, such as permafrost environments with a mixture of rock, air, ice, and water. This method assumes that the four phases are randomly distributed within the subsurface (Somerton, 1992; Luo et al., 1994; Glover, 2010):

$$\begin{array}{}\text{(4)}& \mathit{\rho }={\mathit{\rho }}_{\mathrm{r}}^{f_{\mathrm{r}}}\cdot {\mathit{\rho }}_{\mathrm{i}}^{f_{\mathrm{i}}}\cdot {\mathit{\rho }}_{\mathrm{w}}^{f_{\mathrm{w}}}\cdot {\mathit{\rho }}_{\mathrm{a}}^{f_{\mathrm{a}}}.\end{array}$$

However, the PJI-GM model (Eq. 4) encounters challenges due to the introduction of other unknowns, namely the resistivities of rock, ice, and air, into the system of equations, which can be difficult to determine. Further, inclusion of f_i, f_r, and f_a in Eq. (4) increases the coupling between the system of equations, which may potentially impact the inversion convergence, i.e., referring to the process by which the iterative algorithm reaches a stable solution that adequately fits the observed data. It can therefore result in larger misfits (χ² and RMSEs) between measured and modeled data (Mollaret et al., 2020). PJI-GM has not yet been tested extensively; here, we test PJI-GM on the 22 profiles in Central Asia and assess its suitability for different landforms and substrates.

Following Hilbich et al. (2022), to facilitate comparisons of potential ground ice content across profiles, we define a zone of interest (ZOI) for each profile below the active layer from which a mean ground ice content value is extracted. To define the ZOIs, we typically selected a zone below the active layer that extends horizontally across as much of the profile as possible within the area where frozen conditions are expected. The depth and width of the ZOI were adjusted to each profile's resolution capacity in the relevant depth range to ensure a representative selection (Fig. A2).

3.2.2 Model setup and parameter choice

To determine the best-guess ground ice contents for individual profiles, a systematic analysis was conducted using PJI-GM. This involved iteratively testing 450 different combinations of prescribed subsurface resistivity values (ρ_r, ρ_w, ρ_i, and ρ_a) and start porosities (ϕ_start). The start porosity (initially homogeneous across the profile) is a required initial value for the model and is iteratively adjusted during the inversion. The porosity is directly related to the rock content ($\mathit{\varphi }=\mathrm{1}-f_{\mathrm{r}}$). As our study areas include very diverse landforms and substrates, the measured apparent resistivities exhibit significant variability across profiles, and material-specific properties such as ρ_r, ρ_w, and ρ_i are expected to vary as well (e.g., due to a different ion content in the soil and pore water). To address this variability, ranges of resistivity values representative of all landforms were chosen. For example, fine-grained sediments typically contain higher liquid water contents, leading to lower resistivities. Conversely, rock glaciers often display exceptionally high resistivities due to increased ice and/or air contents in the coarse blocky subsurface matrix. The chosen resistivity ranges correspond here to physically plausible ranges found in the literature (e.g., Telford et al., 1990; Hauck and Kneisel, 2008). Table 3 indicates which parameter values were tested for ρ_r, ρ_w, and ρ_i. The resistivity of air (ρ_a) should theoretically be infinite. However, we found that for fine-grained sediment profiles, ρ_a needed to be lowered significantly (i.e., set to 100 000 Ω m); otherwise the inversion would not converge at all and would result in very high χ² and RMSEs. As a result, we chose to treat ρ_a as a tuning parameter for model convergence, under the constraint that it always remains higher than the resistivities of the other fractions.

Table 3 Values for ρ_i, ρ_r, and ρ_w tested in the PJI-GM loop. Combinations were limited to those meeting the condition ${\mathit{\rho }}_{\mathrm{a}}>{\mathit{\rho }}_{\mathrm{i}}>{\mathit{\rho }}_{\mathrm{r}}$ for physical plausibility. Units are in ohm meters (Ω m).

In addition to the varying resistivity values for the different materials, we tested different regularization parameters for inversions using PJI and chose the most adequate parameter values for each profile individually using the L-curve method (following Mollaret et al., 2020; Pavoni et al., 2023; Wagner et al., 2019; Rücker et al., 2017). In addition to the regularization parameters introduced earlier (α and zWeight), PJI includes a volumetric conservation regularization parameter (β). This was set to the default of 10 000 as we found that this value typically leads to satisfactory results regarding mass conservation. Initial testing in our study revealed minimal variations in χ² and RMSE values around this default value. Similarly, Mollaret et al. (2020) and Pavoni et al. (2023) reported limited variation in χ² and RMSE values for β values near 10 000, supporting its use as a robust default parameter. The regularization parameters applied to each profile are listed in Table 4.

Table 4 Selected PJI model regularization parameters for each profile. Parameters m and n are only used in PJI-AR (Archie's law parameters). All other parameters are needed for both PJI-AR and PJI-GM. ${\mathrm{err}}_{{\mathit{\rho }}_{\mathrm{a}}}$ and err_tt are data errors that were estimated for each profile iteratively to get χ² values closer to 1.

Furthermore, we set the minimal rock content to 0.1 and the maximum rock content to 0.9 for all profiles to allow for the detection of bedrock (high rock content) and excess ice (very low to negligible rock content, e.g., in rock glaciers and ice wedges). This range was not restricted further to assess the model performance in an unbiased way and to see how the model performs if the porosity is unknown, which is often the case for remote regions without boreholes. The PJI model parameters used for each profile can be found in Table A1.

3.3 Clustering approach for tomogram analysis

After running the large number of PJI-GM simulations with the different parameter combinations indicated above, a best-guess representation of the subsurface has to be chosen. A quality check was applied to remove models with large misfits based on the χ² and RMSE thresholds, as outlined in the workflow (Fig. 3). The χ² threshold was chosen following guidelines from the literature, where values between 1 and 5 are considered reliable for most applications (Günther et al., 2006) and values up to 10 are acceptable in specific cases (Audebert et al., 2014; Mollaret et al., 2020). Even after applying these thresholds, the number of PJI-GM model runs with satisfactory fits remained large for many profiles. This is a result of the non-uniqueness of the inversion process, wherein multiple subsurface models can effectively represent the same geophysical data. To assess the range of possible subsurface conditions modeled with PJI-GM, we conducted a hierarchical cluster analysis (analogous to the ensemble inversion approach by Rings and Hauck, 2009). This analysis helps to identify and group similar inversion results based on the mean ground ice content and its standard deviation. The clustering steps can be summarized as follows:

1. Feature extraction. From each remaining tomogram realization after the initial quality check (Fig. 3), key features (mean and standard deviation) of ice content distribution were extracted. This feature extraction step facilitated the comparison of the tomograms and grouping into similar clusters (see Fig. 4 for an example).

2. Hierarchical clustering. Subsequently, to analyze similarities across all tomograms and identify overarching patterns, hierarchical clustering was conducted that grouped all tomograms of each profile into three clusters which represented different types of modeled ice content distributions. We performed hierarchical clustering using Ward’s method, which minimizes the total within-cluster variance during merging and is based on the Euclidean distance (L² norm) between data points (e.g., Ogasawara and Kon, 2021). The clustering analysis was conducted using the scipy.cluster.hierarchy module from the SciPy library (version 1.11.2; Virtanen et al., 2020). We chose the same number of clusters (three) for each profile for easier comparison across profiles. The dendrogram plot (see Fig. 4a) visually represents the hierarchical arrangement of clusters, where each individual tomogram is represented at the bottom. The height of the branches indicates the distance or dissimilarity between clusters, reflecting how much variance increases when merging clusters. Clusters merging at smaller distances are more similar than those merging at larger distances, providing insights into the similarity and structure of the tomograms.

3. Cluster averaging and representative tomograms. Finally, the ice, rock, water and air contents from all tomograms within the same cluster are averaged to obtain a single representative result for the four fractions for each cluster. This results in three representative tomograms, corresponding to the three identified clusters, that show different possibilities of the phase distributions in the subsurface.

Figure 4 Dendrogram and scatterplot illustrating the hierarchical clustering of PJI-GM model outputs for the abra02 rock glacier profile (here, for the 53 remaining tomograms after the quality check). (a) Dendrogram resulting from the hierarchical clustering of the extracted features. The x axis represents the number of each tomogram (from 1 to 53) included in the clustering. (b) Scatterplot of mean ice content and the ice content standard deviation with points colored according to their respective clusters, demonstrating the differences in mean ice content and variability among the clusters.

Figure 5 Example of the PJI-GM clustering results for profile abra02 on the Abramov rock glacier. (a) PJI-GM results (ice, rock, water, air contents) for each cluster, as well as the transformed resistivity (ERT) and P-wave velocity (RST) distributions calculated from these distributions. (b) Virtual borehole plots which, for easier comparison, represent the four fractions of all PJI-GM clusters at one point (x=260 m) with depth.

The clustering approach allows for comparing the variability in ground ice content distribution across multiple model runs. Grouping similar tomograms into clusters representing different subsurface models based on the geophysical data facilitates analyzing the large number of PJI-GM model outputs. Figure 5 exemplifies this with three identified clusters for profile abra02, showcasing their tomographic visualizations (ice, rock, water, and air fractions, along with ERT and RST tomograms). Additionally, Fig. 5b provides a virtual borehole plot at x=260 m, illustrating the depth profiles of these subsurface fractions. To determine the best-guess ground ice content estimate from these PJI-GM clusters, we implemented a multi-step approach. Initially, unrealistic subsurface characteristics were identified and flagged for exclusion. These characteristics, with their corresponding rejection thresholds, are detailed in Table 5.

• Abrupt unrealistic transitions. These are sudden and unrealistic shifts in the vertical distribution of any of the four subsurface components (rock, ice, water, air) or transformed ERT or RST tomograms that lacked a clear justification based on the site conditions.

• Implausible air content. This is air content exceeding reasonable limits for the given soil type and depth. For instance, high air content at depth within fine-grained sediments, where compaction and saturation are expected, would be flagged as unrealistic.

• Implausible water content. Excessively high water content at depths where drainage would likely occur was also considered unrealistic.

• Implausibly low rock content. This occurs at the surface and at the depth of fine-grained sediments.

Figure A3 provides an overview of this workflow, including notes on which cluster of each profile was rejected and the rationale for those rejections.

Table 5 Threshold ranges for cluster rejection based on unrealistic subsurface fractions at different depths. The top 5 m was excluded from this analysis, which in most cases excludes the active layer. CB, coarse blocky (RGs, some MOs); FG, fine-grained (SEDs, TSs).

Figure 6 (a) Correlation plot showing the relationship between the inverted resistivity values obtained from individual inversion of the original ERT (ERT-conv) measurements and the transformed resistivity (ERT-transf) values from the different PJI-GM clusters (mean value of all tomograms of one cluster). The term “transformed” refers to the back calculation of resistivity from the inversion results produced by the PJI-GM model. In this example (profile abra02), Cluster 2 can be excluded as it has the poorest fit. (b) Transformed ERT tomograms of all clusters (C1, C2, C3) versus the tomogram resulting from the individual inversion (Sect. 3.1).

This assessment was further validated by site-specific geomorphological observations and, where available, borehole data. For instance, the presence of landforms indicative of high ice content (e.g., active rock glaciers) or other indicators of ground ice presence (e.g., gelifluction lobes, thermokarst, visible ice outcrops) were considered. While borehole data provided valuable validation at the Kumtor site, such data were not available for the other sites. When multiple plausible clusters remained, we conducted a correlation analysis, comparing the transformed ERT from each PJI-GM cluster with the conventional ERT inversion. The cluster with the highest correlation, indicating the best fit to the measured data, was selected (Fig. 6). This ensured the chosen model best reflected the actual subsurface conditions captured by the ERT measurements. Transformed RST tomograms were not used due to their similarity across clusters. We acknowledge that scarcity of validation data in the region, with borehole data limited to the Kumtor site, necessitates a degree of expert judgment in this selection process.

4 Results

4.1 Permafrost characteristics of the different sites and landforms in Central Asia

We collected 38 ERT measurements from 10 distinct locations to assess permafrost presence and qualitatively estimate ground ice content at each profile, distinguishing between ice-rich and ice-poor conditions. Additionally, 22 of these profiles were co-located with RST measurements, where the PJI method was applied. Our initial examination focused on identifying potential patterns in resistivity across the study sites and landform classes, with the aim of generalizing resistivity signatures by landform. Figure 7 shows the resistivity signature of each profile, grouped by the different landforms. The violin plots represent the full range of resistivity values present within each profile and landform, visualizing both the distribution and density of these values. The width of each violin at any given point reflects the relative frequency of that resistivity value, with wider sections indicating more common values and narrower sections highlighting less frequent occurrences. These plots emphasize the variation in resistivity within different landforms, with distinct patterns emerging for each. For example, rock glaciers (panel b) exhibit higher resistivity values, often associated with ice-rich subsurface conditions, whereas fine-grained sediments (panel c) display lower resistivity ranges, typical of more water-saturated or fine-grained materials. Talus slopes (panel d) and moraines (panel a) show intermediate resistivity distributions, potentially reflecting a mixture of rock, ice, and potentially unfrozen materials. Rock glaciers exhibited the highest resistivities (mean = 700 000 Ω m), followed by moraines (mean = 95 000 Ω m), indicating the presence of the ice in the core of the moraine.

Figure 7 Resistivity distribution of all ERT profiles, grouped by landforms and colored by study site. The horizontal dashed line in each subplot indicates the mean resistivity of all profiles within the landform class. (a) Moraine. (b) Rock glacier. (c) Fine-grained sediment. (d) Talus slope. Each point represents individual resistivity value from the inverted resistivity distribution, capturing the full range of resistivity values obtained across the entire profile. The wide sections of each violin indicate more frequent resistivity values, while narrower sections reflect less common resistivity values. This illustrates the diversity of resistivity measurements within each profile, highlighting variations that occur at different depths and positions along the landform.

The lowest resistivities were recorded for the fine-grained sediment sites (mean = 5500 Ω m). Talus slope resistivities (mean = 25 000 Ω m) fell between fine-grained sediments and moraines. The resistivity distributions across the different landforms show no significant variation between the geographic regions studied. These results suggest that resistivity signatures from the ERT surveys appear to be landform-specific rather than site- or region-specific. However, the resistivity ranges within one profile are generally quite large, as indicated by the large vertical spread of the violin plots in Fig. 7.

Figure 8 Examples of interpreted ERT profiles of different landforms. (a) Rock glacier with high resistivities below an active layer of about 4 m; (b) moraine, with high resistivities possibly pointing to buried glacier ice; (c) rock glacier and fine-grained sediments; (d) talus slope; (e) fine-grained sediments; (f) fine-grained sediments, where a borehole confirmed saturated ground ice conditions in the uppermost layers. The dotted blue lines on the surface indicate the locations of the electrodes.

To further understand permafrost occurrences and conditions across the different landforms in Central Asia, all ERT profiles were independently interpreted with a focus on the likelihood of permafrost presence and potential ground ice content. Layers of high resistivity, typically located beneath a lower resistivity active layer, were indicative of permafrost conditions. The dataset suggests widespread permafrost across all sites and landform types studied. In rock glaciers and moraine surveys, particularly high resistivities were observed, indicating ice-rich conditions (Fig. 8). Notably, exceptionally high resistivities in moraine samples suggest the presence of buried glacier ice in many moraines within the Central Asia study sites, as seen in moraine profiles such as GOL05 (Fig. 8b).

In addition to the more commonly studied permafrost landforms such as rock glaciers and talus slopes, we also conducted measurements on numerous fine-grained, partially vegetated sediment profiles that are prevalent on the high-altitude mountain plateaus of Central Asia. A new borehole was drilled at the Kumtor study site, where profile KUM04 is located, in 2022 (see Fig. 8f), revealing frozen conditions until at least the maximum measurement depth of 32 m, a shallow active layer ranging from approximately 1 to 1.5 m thickness, and saturated ground ice conditions in the upper part of the drill core. Visual inspection of the core revealed fine-grained sediments with a notable presence of interstitial ice, although quantifying the exact ice content was difficult. Based on the sediment structure observed, we estimate a porosity of approximately 20 %–25 %, which would imply a maximum ice content of 25 % under saturated conditions, as is the case in the uppermost few meters below the active layer. At greater depths, while the borehole temperatures remain <0 °C, the sediments appear to be much drier, indicating decreasing ice contents with depth. The low resistivity values (around 5000 Ω m) in these fine-grained sediments suggest that either the ice content is relatively low (likely around 20 %, as indicated by our visual interpretation of the KUM04 drill core) or the sediments retain a significant amount of unfrozen water (or both). Comparable resistivity values (Fig. 7c) were observed in other fine-grained sediment profiles at other sites, implying similar subsurface conditions in those areas.

4.2 PJI-GM clustering results

We used a clustering approach to evaluate PJI-GM and select a best-guess cluster. Figure 9 uses virtual borehole plots to illustrate the variability in subsurface composition (ice, rock, water, and air content) between different clusters for four representative profiles and landforms. Across most profiles, ice and rock content showed the greatest variation between clusters, while water content remained relatively consistent. Air content, consistently low across clusters, appears to be well constrained by the model. However, in some profiles, PJI-GM overestimated surface water content and underestimated rock content, resulting in unrealistic water content values exceeding 50 % (Fig. 9c). This discrepancy was particularly pronounced in fine-grained sediment profiles, where higher surface rock content is expected than what was modeled for some clusters.

Figure 9 Virtual borehole plots of all clusters and mean subsurface fractions of four representative profiles/landforms: (a) SUE05 (MO), (b) abra02 (RG), (c) KUM04 (SED), and (d) abra05 (TS). Red arrows mark examples of unrealistic cluster results.

To quantify ground ice content, we extracted the mean ice content within a representative zone below the active layer (ZOI, Fig. A2) for each profile and PJI-GM cluster (Fig. 10). This is illustrated in Fig. 10, where the selected best-guess PJI-GM cluster for each profile is highlighted with a red outline. Clusters that were rejected based on unrealistic results in any of the tomograms are marked with pink crosses. The figure also shows the mean ice content from the PJI-AR runs as a direct comparison. This will be further discussed in Sect. 4.3. Furthermore, the ice content tomograms of all clusters for the four representative profiles/landforms are shown in Fig. 11. In the following, we summarize some of the main results for each of the four landform classes, focusing on the ice contents.

• Moraines. The PJI-GM results for the three moraine profiles show considerable variability in ice content between individual clusters, highlighting the method's sensitivity to parameter choice. For instance, in the abra10 profile, mean ice content within the clusters ranged from 0 % to 55 % (Fig. 10a). The abra10 and SUE05 profiles appear to feature ice-rich permafrost conditions, as evidenced by the visible presence of ice outcrops within the associated moraine deposits. Conversely, lower ice contents were found in older, more distant moraines with finer-grained sediments, like GOL07b. The estimated mean ground ice content in the sampled moraines is between 15 % and 38 %.

• Rock glaciers. Analysis of the six rock glacier profiles using PJI-GM yielded consistent results across clusters, with all clusters indicating the presence of ground ice (Fig. 10b). Best-guess clusters, selected through a combination of correlation analysis and field observations, yielded ground ice content estimates ranging from 40 % to 60 %, effectively reflecting the high ice contents expected in these landforms. The relatively low variability between clusters and the consistently realistic subsurface tomograms suggest well-constrained conditions for the PJI-GM model in the rock glacier profiles.

• Fine-grained sediments. The PJI-GM model exhibited significant variability in ground ice content estimates between clusters for fine-grained sediment profiles. While most profiles yielded mean ground ice contents within the ZOI ranging from 0 % to 20 %, some exceptions exceeded this range, with profile no. 457 showing the greatest variability (10 % to 70 %). This inconsistency highlights the challenge of characterizing ground ice in these environments using PJI-GM. For example, in profile KUM04, where ice-saturated conditions were observed in the uppermost layers, the majority of PJI-GM runs indicated a complete absence of ground ice (Fig. 11). Despite these challenges, best-guess estimates, informed by correlation analysis, borehole observations, and the exclusion of unrealistic tomograms, suggest a range of mean ground ice contents of 0 % to 23 %. However, systematic parameter estimation for this landform class using the current PJI-GM model remains difficult.

• Talus slopes. The PJI-GM model results for talus slope profiles also exhibit considerable variability. The differences in mean ice content between clusters range from 0 % to as high as 63 %. Some model clusters likely overestimated the ground ice content in talus slopes, leading to unrealistic rock fraction tomograms and their subsequent exclusion from the best-guess selection. The best-guess clusters were selected based on the highest correlation between conventional and transformed ERT data, with mean ground ice contents calculated for talus slopes ranging from 20 % to 40 %.

Figure 10 Mean ice content extracted from the zone of interest for each profile, comparing results from the three PJI-GM clusters with those obtained using Archie’s law. The “best-guess” PJI-GM ground ice content result for each profile is emphasized with a red outline. Clusters that generate unrealistic tomograms for any of the other three subsurface fractions (rock, water, or air) or where field observations demonstrate that the cluster does not accurately represent the subsurface conditions, are marked with pink crosses and are therefore not considered valid representations of the subsurface. (a) Moraine. (b) Rock glacier. (c) Talus slope. (d) Fine-grained sediment.

Figure 11 Comparison of the ground ice content results of the three PJI-GM clusters for a moraine (a), a rock glacier (b), fine-grained sediment (c), and a talus slope (d). For each cluster, the number in the top-left corner indicates how many tomograms (or PJI-GM model runs) correspond to this cluster.

4.3 Comparison of Archie's law (PJI-AR) and the geometric mean model (PJI-GM)

We compared the best-guess PJI-GM cluster results with the more commonly used PJI-AR runs, maintaining consistent regularization parameters, to evaluate the suitability of each model version for the four landform classes. While ground truth data are limited to the Kumtor site, the PJI-GM model consistently produced tomograms with more distinct subsurface structures. In contrast, the PJI-AR model yielded more uniform subsurface fraction distributions across all landforms, as evidenced by lower standard deviations in the tomograms (see also ice content boxplots in Fig. A4). This improvement with the PJI-GM model is likely attributed to its ability to more precisely characterize the subsurface porosity distribution.

Figure 12 Comparison of PJI-AR and best-guess PJI-GM model results for all four fractions (ice, rock, water, air) of the four representative profiles introduced in Fig. 9. The arrows indicate the clear overestimation of water contents and underestimation of rock contents in profiles abra05 (TS) and KUM04 (SED) from the PJI-GM model. The standard deviation of each fraction is noted in the top-left corner of each tomogram.

Compared to PJI-GM, PJI-AR consistently produces lower ground ice content estimates for most rock glacier profiles (Fig. 10b), averaging around 25 %. This seems underestimated given previous findings on rock glacier ice content (e.g., Arenson et al., 2002; Monnier and Kinnard, 2013; Scapozza et al., 2015; Bast et al., 2024). The discrepancy likely arises from PJI-AR's difficulty in distinguishing between ice and rock in high-resistivity settings, which are characteristic of rock glaciers with pore spaces filled with non-conductive fluids. While PJI-AR performs comparably to PJI-GM for profile abra04, PJI-GM generally resolves this ice–rock ambiguity more effectively in all the other RG profiles, generating more realistic ice content estimates and delineating subsurface structures more clearly (Fig. 13). The figure shows that the PJI-AR model yields spatially uniform results for both rock and ice fractions, while the PJI-GM model provides a clearer delineation of the active layer and other structures. This is also quantitatively supported by the higher standard deviations for PJI-GM (see Fig. 12 1a vs. 1b), indicating greater variability in the tomograms. Higher standard deviations suggest the PJI-GM model captures more heterogeneous subsurface features. Notably, both models produce similar results for the other subsurface fractions of rock glaciers.

Figure 13 Comparison of ice and rock contents modeled with PJI-AR and PJI-GM for different rock glaciers. Columns 1 and 3 show the ice and rock content modeled with PJI-AR, while columns 2 and 4 show the same but modeled with PJI-GM.

For talus slopes and moraines, the PJI-AR model generally provides more realistic ice content estimates compared to its performance for rock glaciers, with results often aligning well with those from the PJI-GM model (Fig. 10c and d). However, both models exhibit limitations. In cases where field observations suggest more ice-rich conditions (abra10, SUE05), the PJI-AR model appears to underestimate ground ice contents for TS and MO profiles too. On the other hand, the PJI-GM model frequently overestimates water content in areas characterized by low resistivity, particularly in the near-surface layers of talus slopes, moraines, and fine-grained sediment profiles (Fig. 12f, columns 3 and 4). This overestimation is evident in cases like profile abra05, where PJI-GM predicts water content exceeding 60 % at the surface, contradicting field observations. This tendency to overestimate water content often leads to unrealistic results and subsequent exclusion of PJI-GM clusters, especially in SED profiles. In contrast, PJI-AR does not exhibit this overestimation of water content. Otherwise, PJI-AR appears to be more suitable for characterizing ground ice content in fine-grained sediment profiles, particularly when extensive and time-consuming parameter optimization for PJI-GM is not feasible. This is supported by PJI-AR's accurate capture of low ice content in profiles KUM01, KUM02, and KUM04, where PJI-GM often resulted in tomograms with 0 % ice content or failed to converge.

4.4 Model sensitivity analysis

A sensitivity analysis of the PJI-GM model, using 450 different resistivity combinations (ρ_a, ρ_w, ρ_r, ρ_i) for each profile within the PJI-GM model, revealed that ice content estimates are most sensitive to the choice of ρ_i and ρ_r. Increasing these resistivities generally decreased ice content estimates (Fig. 14), except for most fine-grained sediment profiles (Fig. 14c), where no correlation was found. The analysis further revealed that variations in the parameter values of the water resistivity (ρ_w) did not consistently influence the resulting ice content estimates significantly. Nonetheless, it was observed that lower ρ_w values, ranging from 2 to 20 Ω m, generally produced better model performance across all landforms compared to higher ρ_w values, which less frequently led to model convergence (not shown). Notably, for fine-grained sediment profiles, the ice content consistently approached 0 % whenever ρ_w exceeded 20 Ω m.

Figure 14 Influence of ρ_r and ρ_i on resulting ice content for example profiles of (a) moraine, (b) rock glacier, (c) fine-grained sediment, and (d) talus slope.

Figure 15 Influence of start porosity (ϕ_start) on the modeled volume fractions of ice, water, rock, and air for four profiles that we consider representative of the four different landforms (abra02, abra05, SUE05, and KUM04) for both PJI-AR and PJI-GM model versions. Model runs with χ²>10, indicating a poor model fit, are marked with crosses in the colors corresponding to the profile and are excluded from further analysis.

We investigated the influence of start porosity on modeled ice contents, hypothesizing that higher initial porosities might lead to higher final porosities and, thus, larger ice contents. To test this, we ran the PJI-GM model with varying start porosities (20 %–80 %) while keeping all other parameters constant, using median resistivity values from the best-guess clusters for four example profiles (Table 6). Figure 15 shows the mean values calculated across all the tomograms for the four representative profiles. Contrary to our hypothesis, start porosity does not seem to significantly influence ice, water, rock, and air fractions in most cases. For instance, rock glacier profile abra02 showed only a slight increase in mean ice content (35 % to 40 %) up to a start porosity of 50 % but remained constant thereafter. Similarly, profiles abra05 and SUE05 showed no clear increasing trend. However, the analysis highlighted the PJI-GM sensitivity to parameter changes. Profile KUM04 exhibited a poor model fit (χ²>10) for most start porosities, indicating that even minor parameter adjustments can hinder convergence. This highlights the need for careful parameter selection and interpretation of model results.

Table 6 Resistivity values for ice, rock, water, and air used for the porosity sensitivity analysis in Fig. 15 for all four representative profiles.

Figure 15 also illustrates the impact of the start porosity on the PJI-AR model version, revealing that higher initial porosities result in increased mean ice contents, particularly for the fine-grained sediment profile (KUM04) and the rock glacier profile (abra02), making PJI-AR more sensitive to the chosen start porosity. In contrast to the PJI-GM model, on the other hand, where minor parameter adjustments can induce substantial χ² errors, the PJI-AR model demonstrates greater robustness. Specifically, with the exception of the lowest initial porosity (20 %) for profile SUE05, all other porosities lead to model convergence and acceptable χ² values.

5 Discussion

5.1 Permafrost ground ice contents in the Tien Shan and Pamir of Central Asia

In this study, we modeled permafrost ground ice contents for different landforms in the Tien Shan and Pamir of Kyrgyzstan and Tajikistan using two versions of the PJI model which employ two different petrophysical relations for electrical resistivity: Archie's law (PJI-AR) and the geometric mean model (PJI-GM). The modeled ground ice contents provide valuable insights into typical ground ice contents in various landforms within the mountain permafrost of Central Asia, where such information is currently lacking. Our findings suggest that ground ice contents are primarily influenced by landform type rather than geographic region (e.g., Tien Shan versus Pamir). The mean ground ice contents for the investigated landforms are as follows: ice content for rock glaciers, 38 %–60 %; for moraines, 18 %–40 %; for talus slopes, 20 %–40 %; and for fine-grained sediments, 0 %–20 %.

The results indicate ground ice presence across all investigated landforms, including fine-grained (vegetated) sediment profiles, as confirmed by the borehole core at KUM04 (Fig. 16). This profile, located on a high-altitude plateau at approximately 3540 m a.s.l., resembles surface conditions on the Tibetan Plateau (e.g., Buckel et al., 2021; Gao et al., 2016; You et al., 2017), where thermokarst lakes suggest widespread ice-saturated conditions in various stages of degradation, especially in area 2 of Fig. 16. Our modeling at KUM04 shows saturated conditions with low ice content (20 %), though ice-rich conditions may be more common in thermokarst lake areas, requiring further confirmation due to the lack of RST measurements in this specific area. Upscaling the estimated ground ice contents at KUM04 to section 1 in Fig. 16 (12 km² with a 5 m permafrost layer at 20 % ice content) yields a ground ice volume estimate of 120 000 m³. These ice-saturated conditions in fine-grained sediments contrast conditions in the European Alps, where such sediments typically lack permafrost (Hoelzle, 1994; Kenner et al., 2019). Despite the absence of surface geomorphological indicators, these findings highlight the importance of geophysical methods for detecting and characterizing permafrost, especially on high-altitude plateaus like Central Asia. Such areas host critical infrastructure, including the Kumtor gold mine, where permafrost degradation poses risks of ground instability and thermokarst lake formation.

Figure 16 Map showing one of the regions where ground ice is probably extensive (Gruber, 2012) on a high-altitude plateau with fine-grained, partly vegetated sediments. Photograph (1) shows the location of the borehole; photograph (1a) shows a part of the borehole core with visible ice lenses; photograph (2) shows one of the thermokarst lakes in region 2. Sources of the background hillshade map: Esri, DigitalGlobe, GeoEye, i-cubed, USDA FSA, USGS, AEX, Getmapping, Aerogrid, IGN, IGP, swisstopo, and the GIS User Community.

The mean ground ice contents quantified for rock glaciers (38 %–60 %) in our study region fall within the range reported by other studies, where ground ice contents were either empirically estimated or quantified using the 4PM or PJI approaches in the Andes (e.g., Rangecroft et al., 2015; Jones et al., 2019; Schaffer et al., 2019; Jones et al., 2021a; Mathys et al., 2022; Halla et al., 2021; Hilbich et al., 2022) or in the European Alps (Mollaret et al., 2020; Duvillard et al., 2018; Pavoni et al., 2023). Furthermore, drill cores taken from rock glaciers confirm the potential for high ice contents within these landforms (e.g., Vonder Muhll and Haeberli, 1990; Vonder Mühll and Holub, 1992; Monnier and Kinnard, 2013; Krainer et al., 2015).

For other permafrost landforms, data are scarcer. In fine-grained sediment profiles, we found mean ice content of approximately 20 %, consistent with findings in the Andes (Hilbich et al., 2022) using the 4PM (Hauck et al., 2011). Quantitative data on ice content in talus slopes are also limited (Scapozza et al., 2015), but our estimates (20 %–40 %) align with findings from a Swiss Alps study (Mollaret et al., 2020). Although ice-rich moraines have been extensively studied (e.g., Bolch et al., 2019), quantitative ice content estimates are limited. Our findings (18 %–40 %) align with recent estimates of 40 % in the European Alps (Kunz et al., 2022) but suggest greater variability within and across individual locations and landforms (Fig. 10). These comparisons highlight the variability and uniformity of ground ice content across mountain permafrost landforms and ranges.

Finally, the presence of ground ice in talus slopes, moraines, and fine-grained sediments at all study sites highlights the need to consider permafrost landforms beyond rock glaciers. This is crucial for assessing ground ice content and its hydrological significance in different regions (e.g., Azócar and Brenning, 2010; Jones et al., 2019, 2018). The role of permafrost in mountain hydrology remains poorly understood and is often assumed negligible due to limited data (van Tiel et al., 2024). The estimates presented here are a step towards addressing this knowledge gap, particularly regarding cryosphere–groundwater connectivity. In Central Asia, obtaining more validation data would enhance regional ground ice content estimates, but logistical and financial constraints make extensive borehole drilling impractical. Thus, non-invasive geophysical methods, despite their uncertainties, are critical for advancing permafrost research in this direction.

5.2 Evaluation of PJI-GM for different landforms

The ground ice content estimates presented were modeled using two different versions of the PJI approach. We evaluated the suitability of the PJI geometric mean model, which has not been extensively tested previously, for quantifying ground ice contents across various landforms. While the geometric mean model can generally be applied to all the sampled landforms using the proposed methodology, differences were observed between the landforms and across multiple model runs in terms of the estimated ground ice contents and their spatial distribution within each profile. This is exemplified by the distinct PJI-GM clusters identified, particularly for finer-grained landforms.

The PJI-GM model appears to be most suitable for profiles with a distinct ice-rich layer, as observed in rock glaciers and ice-rich moraines. The examination of ERT tomograms indicated that all rock glaciers display a similar resistivity distribution pattern, characterized by an active layer overlying a thick, high-resistivity layer, which differentiates them from most other profiles where the spatial resistivity distribution is more variable across the longitudinal profile. In these cases, the ambiguity between different PJI-GM cluster results and derived subsurface ground ice contents is minimal. This specifically good performance of the PJI-GM for ice-rich permafrost profiles compared to the PJI-AR model is likely attributable to the fact that the PJI-AR model does not include the resistivity of rock and ice, so the ice–rock ambiguity is solely constrained by the seismic (RST) data.

Figure 17 (a) Ice content boxplots for all rock glaciers with PJI-GM and PJI-AR, extracted from a representative zone (ZOI); (b) P-wave velocity boxplots of each profile from the same ZOI (PJI-AR, PJI-GM, and individual inversion); (c) resistivity boxplots of each profile from the same ZOI (PJI-AR, PJI-GM, and individual inversion). The dashed line in (b) marks the P-wave velocity of ice (3500 m s⁻¹).

The PJI-AR version tends to underestimate ice contents of presumably ice-rich landforms. As discussed before, the only exception, where PJI-AR produces results for a rock glacier site similar to those of PJI-GM, is profile abra04. In this specific profile, a massive ice-rich layer was observed through outcrops in the field. Figure 17b shows that the P-wave velocities (both measured and modeled) for this profile are close to the P-wave velocity of ice (3500 m s⁻¹, marked by the dashed line in Fig. 17) with little spatial variability, which sets this profile apart from other rock glacier profiles where mean P-wave velocities are usually either higher or lower or show a larger spread within the representative zone (ZOI). These results suggest that rock glacier abra04 may contain a core of pure ice from remnant glacier ice, distinguishing it from the other sampled rock glaciers. For this profile, combining Archie's law with Timur's equation likely helps to model the higher ice content, even though electrolytic conduction is improbable. This contrasts other rock glaciers where P-wave velocities are typically higher, closer to those of rock than ice. In this case, PJI-AR has more difficulties in distinguishing between rock and ice. Interestingly, even when seismic velocities exhibit only minor fluctuations around the 3500 m s⁻¹ threshold, as observed in profile abra09 or abra03, PJI-AR results converge to corresponding ground ice contents that are significantly lower than those of PJI-GM. Therefore, we would prioritize using the PJI-GM model for rock glacier profiles, where heterogeneous compositions are likely, especially if no boreholes or other information about the subsurface is available.

For less ice-rich landforms with a lack of a distinct high-resistivity layer and where profiles exhibit more variable and generally lower measured resistivities, PJI-GM produces ambiguous results, with clusters that strongly differ regarding the modeled subsurface conditions. The resulting subsurface ice content models vary here significantly across different clusters (see Fig. 10). For talus slopes and certain moraines analyzed, the fine- to medium-grained materials may facilitate increased subsurface water storage. This higher water content within these landforms may lead to enhanced electrolytic conduction, possibly accounting for the more comparable results between the PJI-GM and PJI-AR models. While the modeled talus slopes and some moraines in this study exhibited relatively fine-grained materials, this may not be representative of the general grain size characteristics of these landforms across different sites. Talus slopes and moraines can display significant variability in grain size, often featuring coarse blocky compositions comparable to rock glaciers, depending on factors like local geology. In such cases, where coarser materials dominate, the PJI-AR model might also underestimate ice content, which is similar to the limitations observed for rock glaciers. This suggests that the PJI-GM model may be more appropriate for these coarser-grained settings. Unfortunately, the available data from this study do not allow for confirmation of this assumption.

The PJI-GM model exhibited more ambiguous performance when applied to fine-grained sediment profiles. For many SED profiles, it was challenging to find parameter combinations that led to successful inversion convergence. Numerous runs either resulted in a lack of convergence, characterized by high RMSEs and χ² errors, as also described in Mollaret et al. (2020), or produced unrealistic tomograms with improbable subsurface fraction distributions. This difficulty in achieving reliable results with PJI-GM for landforms with lower resistivity, such as the SED profiles examined in this study, highlights a potential limitation of PJI-GM. However, selecting lower rock resistivity values (ρ_r), below the minimum threshold of 1000 Ω m used in this study, could potentially help reduce the unrealistic water content estimations observed in some of those profiles. Furthermore, estimating ground ice content in fine-grained sediments at our study sites is inherently challenging due to the absence of clear surface expressions like outcrops or gelifluction features, which are typically absent from flat terrain (Matsuoka, 2001). Given these factors, the PJI-AR model might be generally better suited for analyzing fine-grained sediment profiles, particularly when a priori knowledge of ice/rock content is available. This can also be justified by the higher probability of electrolytic conduction being the dominating electrical conduction process in these fine-grained sediment profiles (e.g., Mele et al., 2014). Furthermore, our findings indicate that PJI-GM has a tendency to overestimate water content in certain talus slope, sediment, and moraine profiles for layers with low resistivities. This leads to unrealistic water content estimates in the tomograms (>50 %–60 %) in layers with low resistivities, as frequently observed in the uppermost layers of our fine-grained sediment profiles (active layer). This overestimation of water content by PJI-GM in low-resistivity layers is likely due to the geometric mean calculation heavily weighting the low resistivity of water (see Eq. 4). Each component's resistivity (ρ_r, ρ_i, ρ_w, ρ_a) is raised to the power of its respective fraction. Due to the significantly lower resistivity of water (ρ_w), the geometric mean calculation tends to overemphasize its contribution in low-resistivity layers, leading to an overestimation of water content.

While constraining maximal water fractions within the PJI-GM model could potentially mitigate this overestimation, our initial attempts revealed a tendency for such constraints to hinder inversion convergence. One potential solution to this issue could be to include a penalty function or to modify the weighting scheme in the geometric mean calculation to reduce the disproportionate influence of low-resistivity phases like water. This approach would involve adjusting the model to more accurately reflect the physical conditions of the subsurface, potentially improving the reliability of the PJI-GM model in low-resistivity environments. Furthermore, accounting for surface conduction, as noted by Mollaret et al. (2020) and Steiner et al. (2021a), for example, might further improve the results.

5.3 Uncertainties in the ground ice content quantification of this study

The uncertainty in the quantification of the ice content presented depends on several factors, which are discussed below:

• i.

  Standard uncertainties in the geophysical data such as the general ERT and RST data quality, resolution capacity, and investigation depth of the surveys, as well as potential inversion artifacts, can all impact the individual inversions, as well as the PJI modeling results. Poor data quality leads to errors in the inversion process, impacting the accuracy of the resulting subsurface models (Hilbich et al., 2009; Hauck et al., 2011; Hilbich et al., 2022; Wagner et al., 2019; Mollaret et al., 2020).

• ii.

  The PJI model relies on numerous parameters that can influence the model output and, consequently, the estimated ground ice content. For example, Mollaret et al. (2020) and Wee et al. (2024) suggest that higher start porosities typically lead to higher ice contents. However, our findings indicate that the initial porosity does not systematically impact the resulting ground ice content, particularly when using the PJI-GM model. Here, increasing, decreasing, and indifferent ice contents were observed with increasing start porosity, with maximum changes in ice content of 10 % (Fig.15). In our study, we found that the start porosity has a larger impact on PJI-AR results than on PJI-GM results (maximum changes of 45 % for the rock glacier profile abra02). Furthermore, incorporating constraints on subsurface fractions, such as limiting the range of rock content (porosity) or water content, could further enhance the model accuracy. This can be particularly useful when detailed a priori information about the subsurface is available. Constraining the rock content can be particularly valuable in cases where the unconstrained model struggles to accurately reproduce observations or when significant ambiguity exists between PJI-GM clusters. For example, PJI-GM overestimates ground ice content in profiles like SUE02 (Cluster 2) and no. 457 (Cluster 3), leading to unrealistic tomograms for the rock fraction. While adding constraints on subsurface fractions could improve model accuracy, this requires detailed a priori knowledge of subsurface conditions, which is often unavailable for remote study sites lacking borehole data. Due to the remote locations studied and the intention to assess the model without bias, we opted to not systematically apply additional constraints to the subsurface fractions.

• iii.

  Equation (2), which links P-wave velocities to the volumetric fractions of ice, rock, air, and water, may not be sufficiently accurate for our study. Our findings indicate that very high P-wave velocities (>5000 m s⁻¹) are not well reproduced in the transformed RST tomograms across most profiles and landforms. Incorporating a more sophisticated petrophysical relationship for the seismic data could further enhance the accuracy of estimating the four subsurface phases, leading to more reliable ground ice content quantification.

• iv.

  The lack of ground truth data presents a significant challenge in validating ground ice content results. Our hierarchical clustering approach for analyzing the PJI-GM model outputs mitigates some uncertainties by showing different possible subsurface ice content distributions, which varied in ambiguity depending on the landform. This method allowed for the exclusion of clusters with unrealistic tomograms in any of the subsurface fractions or transformed ERT or RST tomograms. Additionally, the correlation analysis used to select the final best-guess cluster in most cases coincided with what we would have estimated based on expert knowledge alone. This approach effectively aligns clusters with expert knowledge for most profiles and eliminates implausible results, although it still depends on initial assumptions and the quality of the input data. The variability observed among the plausible clusters (which were not excluded because of unrealistic tomograms in any of the subsurface fractions) can be leveraged to establish uncertainty ranges for the estimated ground ice contents.

6 Conclusions

This study provides a comprehensive electrical and seismic dataset that offers significant insights into permafrost occurrences and ground ice contents in various landforms within the Tien Shan and Pamir regions of Central Asia. By evaluating the applicability of the PJI-GM model for quantifying ground ice contents from geophysical data, we highlighted its strengths and limitations across different landforms.

Our findings underscore the effectiveness of the PJI-GM model for modeling distinct subsurface structures, particularly in ice-rich landforms such as rock glaciers. The PJI-GM model outperforms the formerly used PJI-AR model in these contexts, as it more accurately reflects high ice contents and sharp transitions from low to high resistivity with minimal ambiguity between model runs. This indicates that the PJI-GM model is particularly suited for characterizing ice-rich landforms with layers of high resistivity, making it a more reliable tool for these specific applications. However, PJI-GM encounters challenges with fine-grained sediments, where reduced model convergence and overestimation of water content were observed. In these cases, the PJI-AR model might be more appropriate, especially when extensive parameter tuning is not feasible.

The mean estimates of ground ice content for the various sampled landforms in our study sites in Central Asia can be summarized as follows: rock glaciers exhibit ice contents ranging from 38 %–60 %, moraines from 18 %–40 %, talus slopes from 20 %–40 %, and fine-grained sediments from 0 %–20 %. These results emphasize the necessity of comprehensive assessments across different permafrost landforms to achieve accurate ground ice estimations, beyond just focusing on rock glaciers.

The impact of this study extends beyond the immediate findings. The baseline dataset provided can be instrumental in future monitoring of permafrost dynamics in Central Asia, especially in the context of climate change. Furthermore, this dataset is valuable for modeling studies that aim to assess the hydrological impacts of permafrost thaw and for conducting hazard assessments or adaptation planning for infrastructure built on permafrost.

In conclusion, the PJI-GM model proves to be a valuable tool for specific landforms, particularly ice-rich environments, but requires careful parameter selection for other landforms to ensure accurate results. The insights and baseline data from this study contribute significantly to the understanding of permafrost in Central Asia, laying the groundwork for future research, monitoring efforts, and climate impact assessments. These efforts aim to initiate and partially sustain future endeavors to establish a long-term monitoring network for the essential climate variable (ECV) permafrost, crucial for understanding and mitigating the impacts of permafrost degradation on infrastructure and hydrology in the region.

Appendix A

Figure A1 ERT tomograms (Wenner array) of all profiles, sorted by the study sites. Profile information and statistics are summarized in Table A1.

Figure A2 Individual ERT and RST inversion results and best-guess PJI-GM ice content estimates for all profiles. The zone of interest (ZOI), over which mean values were calculated, is drawn in dotted white rectangles. If a profile spans more than one landform, there are multiple ZOIs, one per landform.

Figure A3 Flowchart depicting the analysis process for selecting the best-guess cluster among the PJI-GM clusters of each profile. The process begins with checking the availability of borehole (BH) information. If BH information is available, the cluster closest to the borehole information is accepted as the best-guess cluster. If BH information is not available, the analysis proceeds without it, rejecting clusters based on various criteria: abrupt unrealistic transitions, implausible air contents, implausible water contents, and implausible rock contents. If multiple clusters remain after these rejections, a correlation analysis between ERT-conv and ERT-transf is conducted to determine the best-guess cluster.

Figure A4 Boxplots showing the distribution of geophysical properties for the zone of interest (ZOI) in all profiles. The first subplot illustrates the ice content extracted from the ZOI. The second subplot presents the P-wave velocities corresponding to the same ZOI. The third subplot shows the resistivity values within the ZOI. Each subplot highlights the variability and range of values for the respective property for all profiles.

Table A1 Summary of individual ERT filtering details for various profiles, including the profile year, landform type, array, spacing, profile length, total number of data points, percentage of data points remaining after filtering, root mean square error (RMSE) percentage, and χ² value for the individual inversions. For the electrode array: W, Wenner; DD, dipole–dipole.

Appendix B

Table B1 Summary of RST filtering metadata, including the profile, year, spacing, length, RMSE, χ², and percentage of reliable first-arrival picks.