Larch (Larix sibirica Ledeb.) and Poplar (Populus laurifolia Ledeb.) Are Moving into the Goby Desert

1. Introduction

Larch (Larix spp.) is the dominant conifer species in Siberia. Larch-dominant communities (Larix sibirica Ledeb., L. gmelinii Rupr. and L. cajanderi Mayr.) are the largest forest formation in Russia, which cover more than 40% of the forested area. The majority of the larch-dominant forests (>95%) are located in Siberia and cover about 70% of the permafrost zone. In severe climatic habitats, they form pure stands.

Due to its high cold-resistance, larch forms treelines in the northern and alpine forest-tundra ecotones. Larch is highly drought-resistant due to its high efficiency of water use [1]. This allows larch to grow at the semi-desert level of precipitation and to populate forest-steppe ecotones. A deep and well-developed root system allows larch to grow on drained soils with poor water and minerals content. However, within permafrost zone, the root system of larch is shallow (up to 30 cm or less).

Larches are anemochoric species with a life span up to 600–700 yrs. (up to more than 1,000 yrs in the north of Eastern Siberia). Larches are the most photophilous, shade-intolerant, and fast-growing conifer species. These species are also well adapted to wildfires and populate post-fire burns with abundant seedlings [2]. In a warmer and drier climate, larch has been considered a potential substitute for precipitation-sensitive Abies sibirica Ledeb., Pinus sibirica Du Tour, and Picea obovata Ledeb.) in southern lowlands of Siberia [3]. Observed climate changes are mostly favorable for larch species within the majority of their range. Thus, since the warming start in the 1970-s, L. sibirica tree growth in general has increased and larch-dominant forest densification is being observed in the Ural and Siberian Mountains [4,5,6]. Larch trees are moving uphill in mountains and northward in the polar tundra [7,8].

However, in the changing hydrothermal regime larch may experience limitation by moisture, followed by growth decrease, which has been documented within larch southern lowland habitats and at high latitudes [5,7,9,10]. Prognostic models suggest that changes in tree ranges and growth will primarily occur between the forest and non-forest ecotones [11]. Thus, larch-dominant forests within forest-steppe ecotone (the Trans-Baikal zone) have experienced thinning and mortality due to increased fire rate [2]. Globally, according to Mamet et al. [12], Larix has experienced northern expansion and southern retraction as species respond to changing climate.

Meanwhile, there are still only a handful of data about the larch growth and vitality at the southern range of its habitat in the changing hydrothermal regime. In this paper, we aimed to analyze the unique larch forest that is growing at the very edge of larch range on the sandy dunes at the Mongolian border, at the northern boundary part of the Gobi Desert. In fact, that site (we called it “Tes-Hem”) is a larch refugium. Together with larch, Populus laurifolia is also present in that area, although that species is less drought-resistant than larch. P. laurifolia is a typical species in Siberia and is often found within river valleys and, more rarely, in the forest-steppe ecotone.

Tested hypothesis: the changing hydrological regime may stimulate growth of Larix sibirica and Populus laurifolia and facilitate their migration into sandy dunes.

We are seeking the answers to the following questions:

• What is long-term larch growth dynamics in the larch refugium?
• What is the age structure of the larch trees?
• What are migration rates of the larch and poplar into sandy dunes?
• What are the GPP dynamics of the vegetation in that site?
• What is the relationship between tree growth and migration, GPP, on one side, and eco-climate variables and relief features on the other?

2. Materials and Methods

The studies were based on field survey data, Terra/MODIS and high-res satellite time series, UAV data, trees growth index (GI) data, and regression analysis between GI and eco-climate variables. Alongside, we modeled the moisture deficit at the end of this century (2080-2100).

2.1. Study Area

Study sites were located on the northern boundary of the Gobi Desert. Sites included two larch refugia, i.e., “Tes-Hem” (the main site) and “Kara-Haya” (supplementary site) (Figure 1).

Trees both poplar (mean height ~2 m up to 4 m and mean diameter 6.5 cm) and larch (mean height ~13 m up to 18 m and mean diameter 35-38 cm) formed clusters with closure ~0.4 as well as a separate single trees. Old-growth larches also presented in the multi-stem forms. The mean crown closure was about 0.1. Both species are growing on poor organic sandy soils (Figure 2).

On-ground vegetation formed by shrubs and some grasses with ground cover about ~10-15%. It composed by Berberis sibirica Pall., Cotoneaster melanocarpus Fisch., Caragana bungei Ledeb., Artemisia spp., Agropyron spp., Oxytropis tragacanthoides Fisch., Hedysarum fruticosum Pall., Thesium tuvinense Krasnob.. Rare small moss communities (Rhytidium rugosum (Hedw.) Kindb. + Polytrichum sp.) located in the larches shadow.

2.2. Ground Survey Data

Fieldworks were conducted in 2025 y. Temporary test sites (N = 12) with radius of 9.8 m (S ~0.03 ha) were established within the elevation range 1100–1240 m a.s.l. We describe relief features (elevation a.s.l., aspect and slope steepness). We obtained geobotanical and soil type data. We measured tree height and DBH (diameter at breast height, 1.3 m), described tree physiognomy and stand closure, regeneration density and vitality. For dendrochronological analysis, trees were randomly selected around the area of approximately 1.0 ha. Samples (cores) were taken at DBH height or at the root collar level with an increment borer. In total, 77 larch and 17 poplar trees were sampled.

2.3. Eco-Climate Variables

Eco-climate variables were obtained from the ERA5-Land database and from the nearest weather station Erzin (WMO index #36307, 50°16′N, 95°10′E, H = 1102 m, distance to test sites is 10 - 20 km) using the AISORI online database [13]. Air and soil droughts were estimated based on the SPEI (Standardized Precipitation Evapotranspiration Index) index [14] and Palmer’s scPDSI (Self-Calibrated Palmer Drought Severity Index) [15] and soil moisture content.

The SPEI is a proxy for atmospheric drought [16]. The increase in the SPEI indicates a decrease in atmospheric drought, and vice versa. The SPEI was calculated based on precipitation and evaporation data extracted from the ERA5-Land database using the SPEI library (v.1.8.1) [16,17] with scripting in RStudio v.2025.09.1 build 375 [18].

The scPDSI was calculated in R-Studio using the R v.4.5.1 programming language [19] and scPDSI library (v0.1.3) [20]. The input parameters were total monthly precipitation and potential evapotranspiration which were extracted from the ERA5-Land database [21]. Negative values of scPDSI indicate drought conditions, and positive ones indicate wet conditions [15].

We estimated the total water content below ground by using equivalent water thickness anomalies (EWTA) obtained from GRACE (Gravity Recovery And Climate Experiment). The monthly EWTA data were downloaded using NASA GRACE(-FO) Data Analysis Tool [22] for the period of 2002–2024 with the spatial resolution of one degree. EWTA is a measure that quantifies changes in the total amount of water stored on land, including soil moisture, snow, and groundwater relative to the 2004–2009 time-mean baseline [23]. These gravimetric data represent the vertical thickness of a layer of water that would result from the measured change in terrestrial water mass [24]. The EWTA provides a valuable tool for monitoring drought [25,26], the water cycle [27] and its impact on vegetation [28,29,30]. Hereafter, we shall use the term “total water content” (TWC) instead of EWTA.

2.4. Moisture Deficit Prognosis

As a proxy of the moisture deficit, we used the annual difference (Δ) between the total annual precipitation (PRE) and evaporation (E). We applied monthly climate data extracted from CMIP6, based on the CNRM-CM6-1-HR model [31] with a spatial resolution of ~0.5 º × 0.5 º to estimate PRE and E. This model was chosen because it values are close to the average of all models presented in the IPCC WGI Interactive Atlas tool [32]. We analyzed data corresponding to climatic scenarios SSP4.5, SSP7.0 and SSP8.5 [33,34]. We fitted the forecast data (PRE and E) to the historical data extracted from the ERA5-Land database using the simple shift method, i.e., the forecasted monthly values were shifted by an amount equal to the average difference between the forecast and actual data for 2015–2024 (the overlap period between the forecast and historical data). We calculated and compared Δ between 2000–2024 and 2081–2100.

2.5. High Resolution Remote Sensing Data

We used the UAV (Unmanned Aerial Vehicle) DJI Mavic 3 Multispectral to obtain very-high resolution scenes during the fieldworks conducted in July 2025. The CMOS RGB + G/R/RE/NIR cameras are installed on this UAV. The scenes were taken at an altitude of about 80 m above the ground surface with about 50% overlap. In total 1300 scenes were obtained. We preprocessed the UAV data using the Agisoft Metashape software [35].

We also analyzed hi-resolution data obtained from freely available web-services [36] with a resolution of 2.0–3.7 m for 2006 and 2024 (QuickBird, WorldView-3, Pleiades Neo) to estimate vegetation cover dynamics. Changes in the vegetation cover were determined using ESRI ArcGIS software based on the expert interpretation and semi-automatic classification procedures.

2.6. Dendrochronology Analysis

The total sample size included cores from 77 and 17 larch and poplar trees, respectively. Tree ring widths were measured with the precision of 0.01 mm. Missing rings were detected by the cross-dating method. The quality of cross-dating was checked by the COFECHA 6.02 software [37,38]. Age trends were eliminated by negative exponential or negative linear regression methods (using the ARSTAN 6.02 software [39,40]. Standard and residual tree ring chronologies were generated. Standard chronologies are biweight robust mean values of indexed raw chronologies. Residual chronologies were generated from standard chronologies by eliminating the autoregressive component (using AR modeling). Tree-ring chronologies (in mm) are a proxy of the annual radial increment. We used a unitless growth index (GI) as the metric of tree growth increment. GI is a normalized tree-ring chronology with an average of 1.0 and relatively constant dispersion. We used the residual chronologies in the dendroclimatic analysis [41].

2.7. Vegetation Productivity Data

The gross primary productivity (GPP) of vegetation cover was estimated using the MODIS MOD17A3HGF product (2001–2025) [42]. These data represent raster composites of the annual GPP values (kgC/m²) with a spatial resolution of 500 meters [43]. Time series data were obtained from the EarthData geoportal [44]. The original GPP data were converted into multi-band images. The average values were determined for specific areas within the 250-meter buffered zones around the test sites using the exactextract algorithm [45]. This algorithm is implemented as a Python library for extracting and summarizing the values in portions of a raster dataset [46].

The GPP trends were calculated for each pixel based on the Theil-Sen algorithm using the pyMannKendall library (v.1.4.2) [47]. It is a non-parametric method based on the Mann-Kendall Test used to analyze time series data for consistently increasing or decreasing trends [48]. Spatial changes in the GPP values between 2022-2024 and 2001-2003 were estimated using Mann-Whitney U-test implemented in the SciPy Python library [49].

2.8. Statistical Analysis

We used StatSoft Statistica v.10 [50], R 4.5.1 [51] along with RStudio v.2025 [52], and Microsoft Excel for statistical analysis. The geospatial analysis was performed using ESRI ArcGIS software and Python. The GPP trends were calculated based on the Theil-Sen algorithm [48] using Python. The changes in GPP were estimated using the Mann-Whitney U-test. Trends in the total water content (EWTA) and the growth index were assessed based on simple linear regression using standard methods. The difference between forecasted and historical climate was estimated using the t-test.

3. Results

3.1. Eco-Climate Variables Dynamics

Long-term improvement of the moisture regime (precipitation increase together with the decrease in air and soil droughts since the 1950s) shifted to the worsening of moisture supply since c. 1990 (Figure 3 a-c). Gravimetric measurements of the total water content indicated its long-term decrease in the 21st century (Figure 3e). However, the moisture regime improved during the last decade (Figure 3 a-c, e). Both summer and annual air temperature have been increasing since c. 1990 (Figure 3d).

The wind speed, an essential variable that influences tree vigor, reached 3.2 m s^–1 during winter with about 0.6±0.1 m s^–1 during the growing period (May-September). The prevailing winds blow from the west and northwest (Figure 4).

3.2. The Chronology of Tree Growth Index

Growth chronologies of L. sibirica for Tes-Hem and Kara-Haya sites were merged because both chronologies were strongly correlated (r = 0.68) (Figure 5a).

The larch population consists of two cohorts: old-growth trees (age up to 500+y) and young generation (age up to 60 y). The formation of the young cohort corresponded to the period of air and soil drought decrease in 1980s (Figure 3 b,c; Figure 5c).

Larch growth chronologies indicated warming in the 16th century followed by the Little Ice Age in the 17th and subsequent periods of warming and cooling (Figure 5a). An interesting phenomenon observed in the year 1710, namely a sharp decrease in larch growth. The only plausible explanation for this is damage to the trees caused by a lightning strike.

An air and soil drought decrease in 1980-s –1995-s stimulated larch growth with a following depression until c. 2003, followed by GI increase during the last decades (Figure 5b). The growth of the young larches was similar to that of old-growth trees (Figure 5c).

The poplar tree chronology is more stochastic in comparison with the larch chronology. Similar to larch, poplar growth decreased in 1990s and then increased in the 21st century (Figure 6).

3.3. Trees Growth Dependence on Climate Variables

Growth of both old-growth and young larch tree cohorts is controlled by the moisture variables – precipitation, air and soil droughts. Notably, there is a higher dependence of the young larch cohort in comparison with old-growth trees (r = 0.45…0.75 versus 0.30…0.49) (Figure 7 b-d). This is explained by the differences in the roots habitat zone which is evidently deeper in the case of old-growth trees. Air temperature impairs young larch growth whereas its influence on the old-growth trees was insignificant (Figure 7a).

Poplar trees dependence on air temperature and moisture variables is similar or stronger (in the case of soil drought) in comparison with young larch trees (Figure 7).

Summer wind (July) impaired larch growth presumably by wind-caused water loss via increasing evapotranspiration (Figure 8).

3.4. Trees Proogation into Desert

Surprisingly, the lesser drought-resistance poplar trees migrated into the desert further than larch trees. The average shifts of poplars and larch boundary into desert were 35±5 m (max = 70 m) and 280±10 m (max = 470 m) since 1985 to 2024. The year 1985 is the proxy of the larch and poplar young cohort trees establishment (Figure 9).

The mean migration speed of larch and poplar were 0.9 m/y and 7.0 m/y

3.5. GPP Trends

Positive GPP trends cover about 76% of the sandy dunes area. The rate of GPP increase was ~19 kg/ha/y. Mean GPP in the last five years (2020-2024) increased by 250 kgC/ha in comparison with 2002–2006 (i.e., at the beginning of observations). Currently, mean vegetation GPP is 2080 kgC/ha with max values around 6000 kgC/ha (Figure 10).

3.6. The Moisture Deficit Projections

The moisture deficit (MD) is an important hydrological parameter that indicates the moisture regime of vegetation. MD was calculated as the difference (∆) between precipitation and potential evapotranspiration (PET). We estimated MD up to the end of the 21st century (2081–2100) for the climate scenarios SSP2-4.5, SSP3-7.0, and SSP5-8.5. The scenario SSP1-2.5 was not considered since it is not realistic now.

The increase of MD is predicted in the beginning and middle of the growth period (May–July) according to all climate scenarios. However, changes in the annual MD are insignificant (Figure 11).

4. Discussion

For forest ecologists, the most significant changes in tree growth, vigor and range expected to be within the forest - non forest ecotones. In this work we considered the larch growth within unique ecotone in which larch refugia are in contact with Mongolian desert or, in broader view, with the northern boundary of the Great Goby desert. The age of alive larch in refugia exceeds 500 years which suggests that this “last larch frontier” has existed throughout the Holocene. Refugia is located at the edge of larch range along the geographical meridian. Trees survival in refugia is facilitated by site relief features. Namely, trees grow on sheltered downwind north-west and south-west slopes and within local depressions. These are the zones of the snow accumulation and lower evapotranspiration.

Larch population in the refugia includes old-growth and young tree cohorts (age A = 50+). Young cohort establishment corresponds to the end of 1970th, when warming coincided with local minimum of air and soil drought (Figure 3 a-d). Then, when warming-driven water stress increased in 1990s, both tree growth and regeneration recruitment decreased, followed by an increase in the 21st century (Figure 5). In the 21-s century, the population of L. sibirica increased by +36%.

Larch growth response to a warming climate included several phases. First, since warming began in the 1970s, growth increased. Then, growth stagnation and decrease were observed in the late 1990s and very beginning of the 2000s, associated increased water stress. Finally, in the 21st century larch growth rose against the background of decreasing air and soil drought (Figure 5).

Indeed, larch growth is strongly controlled by the moisture variables – air (SPEI index) and soil (Palmer’s scPDSI index) droughts (Figure 7c, d). Old-growth larch is less sensitive to moisture stress than young larch, which is evidently explained by its deeper root habitat. Young trees are also more sensitive to the negative influence of temperature at the beginning of the growing period (Figure 7a). This influence is mostly indirect, i.e., caused by increased water loss. However, warming has stimulated larch growth in the southern Siberian and northern China mountains [5,9].

It is notable that summer winds impaired larch growth, whereas the influence of winter wind was insignificant (Figure 8). This influence explained by increased water loss trough evapotranspiration. Typically, winter winds impair trees growth due to desiccation and snow abrasion [53]. However, in the study area winter wind speed is low (~3.0 m/s) and its maximal values are decreasing (Figure 4b, c). In addition, larches are protected by the dune’s slope against prevailing winds (Figure 9).

Similarly to larch, relict Scots pine (Pinus sylvestris L.) at its southern edge is also showing increasing growth and regeneration density in the 21st century (Figure 12a) [54,55]. This pine refugia is located northward in the Tuva hollow. Both species are growing on the sandy dunes. Pine is also a drought-resistant species, although its resistance to moisture stress is lower (Figure 12 b, c). Therefore, the pine southern boundary is located at c. 100 km northward. Unlike larch, relict pines have experienced an increasing wild fire influence that has damaged mature trees and, more importantly, burned regeneration. So, the existence of this unique refugia depends on the combination of available moisture and burning rate.

Opposite to reported shrinkage of larch range at its southern edge (e.g., in the Trans-Baikal zone and in the West Sayan Mountains, [2,5]), on the boundary with the Mongolian desert larch is expanding its range. However, the above-mentioned range shrinkage is caused by an increased fire rate, whereas the desert is a no-fire area. Meanwhile, larch forest dieback has been described in the Mongolian Mountains [56]. Earlier studies also revealed larch and Scots pine growth decrease in the Asian arid zone in the beginning of the 21st century [57,58]. However, our longer chronology indicated that both species increased growth during last two decades (Figure 12a).

Within the larch refugium a single old-growth poplar trees (A > 200 y) were found (Figure 2c). The young poplar trees (A = 40+ y) are located southward of the refugia (Figure 9). Similar to young larch trees, they started to establish since warming began in the 1970s. Growth of poplar trees is more strongly controlled by air and soil droughts in comparison with larch trees (Figure 7 c, d). In the 21st century, the population of P. laurifolia increased by +46%.

Both drought-resistant larch and less resistant poplar were slowly migrating into the adjacent desert. Unexpectedly, poplar trees migrated ahead of larch (mean shift and speed are 280 m and 7 m/y vs 35 m and 0.9 m/y for larch). The mean speed of the larch trees migration was about 0.9 m/y, which is about twice as high as the rate of larch and other conifer migration into alpine tundra [53]. Poplar migration speed is about one order higher (7 m/y). Populus propagation facilitated by its high seed productivity (> 1 million seeds for a mature tree). Light seeds with “parachutes” spread over hundreds of meters. In addition, individuals can also reproduce through the roots by suckers. Because of that, poplar regeneration reached 900 seedlings/ha, whereas larch regeneration is sporadic. Therefore, the new niches created by the changing moisture regime are colonized by P. laurifolia ahead of larch.

Trees growth increase coincided with GPP increase of vegetation (composed by small bushes) in adjacent sands. That represents a “greening” of the Gobi Desert at its northern boundary. Similar phenomenon has been observed in the nearby Tuva hollow, where steppe communities are increasing their GPP [54,55]. It is worth noting that described increase in population of trees and shrubs has facilitated sands conservation and suppressed desert northward migration.

In general, we observed a unique event of southward migration of trees at the edge of their range, whereas in number of publications describe tree species range shrinkage at their southern boundary (e.g., [11,59]). For example, P. tremuloides Michx. has experienced negative climate-driven changes in the North American forests [60]. Larix gmelinii decline and mortality have been described in the Trans-Baikal zone [12,61], Betula sp. in Mongolia and the Trans-Baikal zone [62,63]. Conifer mortality at their southern range has been described in Europe and America [64,65,66]. However, predicted increase of moisture deficit and worsening of the hydrological regime (with horizon to 2080–2100) suggests a decrease of trees ranges in the studied zone. Nevertheless, the changing climate is currently favorable for the growth of trees and bushes at the northern boundary of the Great Gobi Desert.

5. Conclusions

At the southern edge of larch and poplar ranges we observed a unique southward migration of trees, while published studies revealed the shrinkage of tree species at their southern ranges throughout the boreal biome. Moreover, trees are migrating from refugia into sandy dunes. Notably, P. laurifolia migration rate is about one order of magnitude faster than that of L. sibirica. Alongside that, GPP of vegetation cover in the desert is increasing. Currently, climate warming is favorable for the growth of trees and bushes at the northern boundary of the Great Gobi Desert.