Study area and data
Study area and validation site
The Qinghai-Tibet Plateau, located between 26°00’ to 39°47’ N and 73°19’ to 104°47’ E, is a vast and geographically diverse region. Spanning ~2800 km from east to west and ranging in width from 300 to 1500 km from north to south, it encompasses a total area of around 2.5 million square kilometers. With an average altitude exceeding 4000 m and varying between 3000 and 5000 m, the plateau is the cradle of numerous significant rivers that flow into East Asia, Southeast Asia, and South Asia. This high-altitude terrain is blessed with plentiful solar and geothermal energy resources, as well as extensive permafrost areas. The landscape is predominantly characterized by natural grasslands, which are diverse and can be classified into eight distinct types: alpine vegetation, alpine shrub meadows, subalpine sparse forest shrub meadows, alpine steppes, alpine mountain deserts, alpine swamps, mountain shrubs, and sparse and mat-like grasslands. Among these, the alpine vegetation stands out as the most extensive and significant grassland type27.
The climate across the plateau is characterized by a sharp gradient in temperature, with the annual average temperature dropping from 20 °C in the southeast to below −6 °C in the northwest. This temperature variation is influenced by the obstruction of warm and humid oceanic air currents by a series of mountain ranges. Similarly, the annual precipitation shows a significant decrease from 2000 mm to less than 50 mm, a pattern also shaped by the region’s topography. On the northern slope of the Himalayas, the annual precipitation is less than 600 mm, while the southern slope experiences a subtropical and tropical north edge mountain forest climate, characterized by an average temperature of 18–25 °C during the hottest month and an annual precipitation ranging from 1000 to 4000 mm. In the southern and western parts of the Kunlun Mountains, the climate shifts to a cold semi-desert and desert climate, with the warmest month’s average temperature hovering between 4–6 °C and an annual precipitation of 20–100 mm. Despite these varied climatic conditions, the plateau benefits from abundant sunshine, receiving an annual solar radiation total of 140–180 kcal/cm² and an annual sunshine duration of 2500–3200 h (China Digital Science and Technology Museum, https://www.cdstm.cn/) .
Overall, the Qinghai-Tibet Plateau is a region of remarkable environmental diversity and ecological significance, shaped by its unique geographical position and complex topography. Understanding its climate and the distribution of its vegetation types is crucial for ecological conservation, sustainable resource management, and the assessment of its role in the global carbon cycle.
Validation data for this study is derived from the eddy covariance measurements at the A’rou Grassland Station, as illustrated in Fig. 14. Nestled on the southern banks of the Bajiao River—a tributary of the Hei River—the A’rou Grassland Station occupies a relatively level and expansive riverplain highland. The station is precisely located at coordinates 100°27’ E and 38°02’ N, with an elevation of 3032 me, where the predominant landscape is that of alpine vegetation.
The comprehensive atmospheric observation program at the station comprises a suite of measurements, including wind temperature and humidity gradients at 2 m and 10 m above the ground, barometric pressure monitoring, and precipitation assessment. The atmospheric, radiation, and soil observations are meticulously recorded at a frequency of 0.2 Hz, with data points averaged over 10-min intervals to ensure the reliability and accuracy of the collected information. To capture more dynamic and transient atmospheric processes, eddy correlation measurements are taken at a higher frequency of 10 Hz. The eddy covariance data utilized in the current research spans the period from 2013 to 2014, providing a rich dataset that enables a thorough analysis of the carbon and energy exchanges within the alpine vegetation ecosystem of the Qinghai-Tibet Plateau.
Data
(1) Meteorological data and processing
The meteorological data that serves as input for the model encompasses a range of daily parameters, including average temperature, precipitation, net radiation, and CO2 concentration. These daily averages for temperature, net radiation, and rainfall are derived from hourly observations, providing a detailed temporal resolution for the model’s analysis. The dataset for daily average temperature spans from January 1, 2008, to December 31, 2014, offering a comprehensive overview of the temperature trends across the Qinghai-Tibet Plateau. The plateau is characterized by generally low daily average temperatures, with the maximum temperatures ranging from 25.0 to 26.0 °C and the extreme minimum temperatures falling below −36.0 °C. In the northern parts of Tibet, the average temperature is notably cooler, with a mean of 8.0 °C. The temperature patterns reveal that the lowest daily average temperatures occur in late January and early February, while the peak temperatures are typically observed in July, August, and September. The daily rainfall data, also covering the period from January 1, 2008, to December 31, 2014, highlights the spatial and temporal variability of precipitation across the Qinghai-Tibet Plateau. Rainfall distribution is notably uneven, with a significant concentration of precipitation in the southeastern region, particularly during the summer months of July, August, and September. The model also incorporates daily radiation data for the same period, from January 1, 2008, to December 31, 2014. This inclusion of radiation data is essential for understanding the energy balance within the ecosystem and for accurately modeling the interactions between the atmospheric and terrestrial components of the carbon cycle. Overall, the integration of these detailed meteorological datasets into the model provides a robust framework for analyzing and predicting the complex dynamics of carbon fluxes on the Qinghai-Tibet Plateau.
To improve the simulation accuracy of the model, it is not only necessary to enhance the functionality of the model itself but also to make full use of existing observational data. The transition of the LPJ model from point-scale to regional-scale applications demands a substantial dataset for input. Generally, the regional input data are derived and upscaled from point observations through various techniques, such as Kriging interpolation, to ensure spatial representativeness. In addressing gaps in the dataset, a suite of strategies is employed. For instance, data points with occasional missing values are often filled by substituting with the mean of adjacent measurements. For instances where more extensive data gaps are present, linear interpolation is applied to estimate the missing values. In cases where there are more than 10 consecutive missing values, the approach involves borrowing values from corresponding time frames of adjacent dates. To generate gridded data at a finer resolution, daily average temperature, net radiation, and rainfall data are subjected to geographic interpolation to fit a 5 km × 5 km grid pattern, using latitude and longitude as the projection basis. This process facilitates a more nuanced and localized representation of the input variables, which is critical for the model’s performance and accuracy.
As for the CO2 concentration data, the LPJ model requires annual totals of CO2 concentration for the period spanning from 2008 to 2014. This essential data is sourced from the Waliguan station located in Qinghai province. The monthly CO2 concentration data collected at the Waliguan station are meticulously processed and aggregated into annual totals, which are then formatted and input into the LPJ model according to the required specifications (https://gaw.kishou.go.jp/).
By meticulously integrating and processing these diverse datasets, the LPJ model is poised to deliver simulations that are both robust and reflective of the complex dynamics at play within the regional ecosystem.
Vegetation remote sensing data
The LAI stands as a pivotal structural parameter that facilitates the quantitative analysis of energy exchange dynamics within vegetative ecosystems. It is precisely defined as half the total one-sided leaf area per unit of ground surface area. In this study, the LAI dataset is sourced from a remote sensing product crafted by Beijing Normal University, which utilizes MODIS (Moderate Resolution Imaging Spectroradiometer) data. The original data boasts a resolution of 1 km, providing a detailed depiction of the LAI across the landscape (https://muses.bnu.edu.cn/). To align with the requirements of this study, the spatial resolution of the LAI data is adeptly adjusted to a 5 km × 5 km grid (Fig. 13). This resizing enhances the compatibility of the data with the model’s input specifications while maintaining a balance between the level of detail and computational efficiency. Subsequently, the data corresponding to the Qinghai-Tibet Plateau is meticulously extracted and transformed into a text format. This conversion is instrumental in streamlining the data integration process into the model, ensuring a seamless and accurate representation of the region’s LAI for the study’s analytical purposes.
The left panel shows the LAI data in 2008, the right panel is that of in 2014.
Land cover and soil data
The land use/cover dataset utilized in this study is derived from high-resolution 30 m LANDSAT satellite imagery (https://data-starcloud.pcl.ac.cn/iearthdata/), refined by the expertise of Tsinghua University (Fig. 14). To harmonize with the study’s requirements, the original data is resampled to a 5 km × 5 km resolution, which optimizes it for the analytical scale of the research. The land use/cover data is meticulously classified into distinct categories to facilitate comprehensive analysis. These categories include bare land, water bodies, forests, grasslands, croplands, and snow/ice. Each of these categories is allocated a proportion (Table 1), with the data indicating that in the Qinghai-Tibet Plateau region, bare land and grasslands constitute the predominant land cover types. Together, they account for 48.7% of the region’s total land area, reflecting the plateau’s characteristic landscape. Furthermore, the aggregated vegetation cover, encompassing forests, grasslands, and croplands, represents 23.81% of the area. This statistic underscores the significant role of vegetation in the ecological functions and environmental health of the Qinghai-Tibet Plateau.
The red pentagram represents the position of A’rou Flux site.
The LPJ model necessitates detailed datasets that capture the intricacies of soil properties, including texture and organic matter content. Soil texture is pivotal as it influences numerous soil characteristics and is categorized into various types, each with its own set of parameters for the model. The LPJ-DGVM takes a comprehensive approach by classifying soil textures into 13 distinct types. This classification is grounded in the global soil texture map at a 1° resolution28. The recognized soil texture types include Heavy-clay, silty-clay, clay, silty-clay loam, clay loam, silt, silt loam, sandy clay, loam, sandy clay loam, sandy loam, loamy sand, and sand. These classifications enable the LPJ-DGVM to accurately simulate soil behavior and its interaction with vegetation across diverse landscapes.
In the context of the Qinghai-Tibet Plateau, the distribution map of soil organic matter content reveals a predominantly narrow range, with most areas exhibiting soil organic matter content between 1% and 2%. This insight is vital for calibrating the LPJ model parameters specific to soil characteristics within the study area. The accurate representation of soil properties is crucial for the model’s ability to simulate key ecological processes such as plant growth, carbon sequestration, nutrient cycling, and overall ecosystem dynamics.
Validation data
The eddy covariance data collected at the A’rou Station spans the years 2013 and 2014, providing valuable insights into the atmospheric fluxes during this period. A 2-h gap in the data on July 14th was encountered, and this interval has been deftly filled using linear interpolation to maintain the integrity of the dataset. It is important to note that certain correction algorithms for eddy covariance data may not be executed during the digitization phase, underscoring the necessity for subsequent data post-processing. The standard workflow for processing eddy covariance data is outlined in the following section. This process is designed to address and rectify various data anomalies that may arise from environmental interference. Outliers, which can significantly skew variance and covariance values, often result from environmental factors such as rain, snow, dust particles, or sudden power interruptions affecting the sensor’s acoustic path.
The protocol for outlier removal involves several steps. Initially, it includes filtering out observations recorded during abnormal data acquisition conditions, as indicated by diagnostic values from the data logger. This is followed by the elimination of observations that fall outside the theoretical maximum ranges for parameters like wind speed, CO2 density, and water vapor density. Subsequently, outliers are identified using the Δx ≥ n σΔ formula, where x denotes parameters such as wind speed, CO2, and water vapor density. In this formula, Δx represents the difference between consecutive observations of the same parameter over time, σΔx is the root mean square of these differences calculated over a specific period, and n is a natural number, commonly set at n = 4. Outliers are then replaced with values derived from linear interpolation to ensure the continuity and reliability of the dataset.
This meticulous processing workflow is crucial for ensuring the accuracy and quality of the eddy covariance data, which is essential for understanding the complex interactions between the atmosphere and the biosphere. The process of coordinate rotation is a critical step in aligning the ultrasonic anemometer’s coordinate system with the natural coordinate system, which is based on the average wind speed components over a specified period. This alignment is essential for fulfilling the assumptions inherent in turbulence correlation analysis. The rotation procedure is as follows: Initially, the w-axis (vertical axis) is kept constant, while the u and v axes (horizontal axes) are adjusted. The u-axis is realigned to coincide with the average wind direction, effectively nullifying the v-axis’s average wind speed to zero. Subsequently, with the v-axis held steady, the u and w axes undergo a rotation that results in the w-axis recording an average wind speed of zero. Momentum flux measurements are particularly sensitive to these coordinate rotations, as they directly influence the accuracy of the flux estimates.
Due to a time lag in the measurement of gas density by the LI-7500 anemometer compared to the ultrasonic temperature readings, this discrepancy is addressed through a lag correction process. This correction ensures that the gas density measurements are synchronized with the temperature data, which is vital for accurate flux calculations.
The Web Processing Library (WPL) algorithm is employed to mitigate the effects of air density variations on the observed trace gas fluxes29. After implementing the WPL algorithm, 30-min data on latent heat and CO2 fluxes are derived. Negative values from the 30-min latent heat data are excluded, and the daily latent heat flux is computed through accumulation. This latent heat flux is subsequently transformed into estimates of evapotranspiration. For the 30-min CO2 flux data, a U* correction is initially applied, with a threshold value set at 0.1 m/s30. Following this, the daily CO2 flux data is aggregated. The entire turbulence data processing workflow is executed within the Edire software environment, ensuring methodological consistency and data integrity.
The processed turbulence data for the year 2013–2014 is presented in Fig. 15, showcasing the comprehensive analysis and the meticulous attention to detail that characterizes this study’s approach to understanding the atmospheric exchanges within the Qinghai-Tibet Plateau region.
The above panel shows the field observation data in 2013, the below panel is that of in 2014.
Overall technical process
By analyzing the distributional traits of the plateau’s vegetation and the original 10 plant functional types of the LPJ model, parameterization of the new functional types tailored to the unique ecological conditions of the Qinghai-Tibet Plateau is carried out to develop an LPJ-DGVM model suited for this high-altitude region. Driven by the latest time-series remote sensing data set, the model is further used to simulate and analyze the impact of climate change on the carbon and water fluxes (GPP, NEP) of the Plateau vegetation. The detailed technical workflow is shown in Fig. 16.
Including four parts: Input Data, Data assimilation and model improvement, Validaton analysis and Impact assessment.
The LPJ-DGVM model, an evolution of the BOIME model series, incorporates dynamic vegetation processes to capture the immediate responses of ecosystems to climate change29. The model’s primary inputs consist of atmospheric temperature, cloud radiation data, precipitation, and other meteorological variables, as well as soil texture and atmospheric CO2 content.
The output includes GPP, NPP, NEP, evaporation, transpiration, and the distribution of carbon and nitrogen in vegetation, litter, and soil organic matter. The LPJ-DGVM simulates 10 types of vegetation, including eight woody species and two herbaceous species. The time scales of the various modules in the LPJ-DGVM model are not consistent, with processes such as photosynthesis, water balance, leaf photosynthetically active radiation (PAR) ratio, sunshine duration, potential evapotranspiration, and soil temperature being on a daily or monthly scale. In contrast, soil texture, CO2 content, maintenance respiration, and photosynthate allocation are modeled on an annual scale. Therefore, the model’s temporal scale can be further improved through data-driven approaches.
LPJ-DGVM model and localization
The simulation of vegetation dynamics is a multifaceted endeavor that encompasses a range of physiological and ecological processes, which can be broadly categorized into four interconnected parts: the carbon cycle, water cycle, energy flow, and vegetation dynamics themselves. To enhance the precision of the model’s simulations, it is imperative to integrate these four dimensions holistically. The LPJ model considers the cycling of materials and the dynamics of vegetation in its framework. It achieves the simulation of populations on each simulation unit by meticulously modeling the physiological processes of average individuals on each unit. These physiological processes can be succinctly distilled into three principal modules: The Plant Photosynthetic Productivity Module, the Vegetation Hydrological Balance Processes Module, and Other Physiological and Ecological Processes Modules. The following sections provide a detailed introduction to the modules that have been refined and improved for the LPJ model, ensuring that the simulations are more accurate and reflective of the complex interactions within the ecosystem.
Vegetation photosynthetic productivity calculation and improvement of spatiotemporal resolution
The environmental factors that primarily affect photosynthesis and vegetation respiration include light, temperature, CO2 concentration, and water. The main purpose of the photosynthetic productivity module is to estimate the NPP of vegetation, which centrally involves simulating plant photosynthesis and plant respiration.
Photosynthesis models include empirical models that are based on the statistical relationship between vegetation productivity and climatic factors, process-based models that are derived from biological and physiological processes, and remote-sensing-driven models that estimate vegetation productivity through light energy utilization by using remote sensing data serving as the primary input. A typical process-based photosynthesis model is the Farquhar model31. The model posits that photosynthesis is mainly limited by PAR, the Rubisco enzyme, and CO2 concentration. It calculates the rate of photosynthesis under each of these constraints and adopts the minimum value as the actual rate of photosynthesis.
Given the capacity of remote sensing to capture surface information, and considering the high temporal resolution of numerous remote sensing satellites, many methods have been developed for estimating vegetation productivity through remote-sensing-driven models. The general form of such models is expressed as NPP = ξ * PAR * FPAR, where ξ represents the light energy conversion rate, which is a function of temperature, moisture, and CO2 concentration. PAR refers to photosynthetically active radiation, and FPAR denotes the absorption fraction of photosynthetically active radiation.
The photosynthesis model in LPJ-DGVM uses a model that seamlessly integrates light energy utilization with Farquhar, allowing the introduction of original daily high-temporal resolution data on environmental factors such as temperature, moisture, and CO2 concentration directly into the photosynthesis model through the Farquhar model. The model’s expression is as follows:
$$\begin{array}{l}{A}_{{nd}}={I}_{d}\left(\frac{\alpha {f}_{{temp}}\left({p}_{i}-{\varGamma }_{* }\right)/\left({p}_{i}+2{\varGamma }_{* }\right)}{\left({p}_{i}-{\varGamma }_{* }\right)/\left({p}_{i}+{k}_{c}\left(1+\frac{{\rm{p}}{{\rm{o}}}_{2}}{{k}_{0}}\right)\right)}\right)\left[\left({p}_{i}-{\varGamma }_{* }\right)/\left({p}_{i}+{k}_{c}\left(1+\frac{{\rm{p}}{{\rm{o}}}_{2}}{{k}_{0}}\right)\right)-\left(2{\rm{\theta }}-1\right)\left(\frac{24a}{h}\right)-2\left(\left({p}_{i}-{\varGamma }_{* }\right)/\right.\right.\\\left.\left.\qquad\qquad\left({p}_{i}+{k}_{c}\left(1+\frac{{\rm{p}}{{\rm{o}}}_{2}}{{k}_{0}}\right)\right)-\left(2{{\rm{\theta }}}^{2}-{\rm{\theta }}\right)\left(\frac{24a}{h}\right)\right){\left[1-\frac{\left(\left({p}_{i}-{\varGamma }_{* }\right)/\left({p}_{i}+{k}_{c}\left(1+\frac{{\rm{p}}{{\rm{o}}}_{2}}{{k}_{0}}\right)\right)-\left(2{\rm{\theta }}-1\right)\left(\frac{24a}{h}\right)\right)}{\left(\left({p}_{i}-{\varGamma }_{* }\right)/\left({p}_{i}+{k}_{c}\left(1+\frac{{\rm{p}}{{\rm{o}}}_{2}}{{k}_{0}}\right)\right)-\left(2{{\rm{\theta }}}^{2}-{\rm{\theta }}\right)\left(\frac{24a}{h}\right)\right)}\right]}^{0.5}\right]{0.088H}_{{est}}^{3}\end{array}$$
(1)
Where, \({A}_{nd}\)(g·cm−2·d−1) represents the daily photosynthesis rate, \({I}_{d}\) represents the total photosynthetically active radiation absorbed per day, and \(\theta\) is a constant. h represents the daylight duration, a is a constant, \(\alpha\) represents the inherent quantum yield of CO2, \({f}_{{temp}}\) represents the temperature inhibition function, Г* represents the CO2 compensation point: \({\varGamma }_{\ast }=\frac{p{o}_{2}}{2\tau }\), \(p{o}_{2}\) represents the oxygen partial pressure, \({p}_{i}\) represents the intercellular CO2 partial pressure: \({p}_{i}=\lambda {p}_{a}\), \({p}_{a}\) represents the atmospheric CO2 partial pressure, \(\lambda\) is a constant, τ、kc and \({k}_{o}\) are dynamic parameters related to temperature32.
These parameters and functions are intricately linked and collectively determine the photosynthetic rate, which is fundamental to the calculation of the ecosystem’s carbon balance and productivity. In models for calculating vegetation NPP, respiration is generally simulated. NPP is calculated as the net difference between Gross Primary Production (GPP) and the total respiration (R), where NPP = GPP − R. Various models simulate respiration differently, with some establishing intricate relationships between total respiration and a suite of environmental factors.
The LPJ-DGVM model adopts a nuanced approach to simulate respiration by distinguishing between maintenance respiration and growth respiration. The LPJ-DGVM adopts a nuanced approach to simulate respiration by differentiating between two primary types: maintenance respiration and growth respiration. The basic idea is that the maintenance respiration of various plant tissues is influenced by phenological status, tissue temperature, and nitrogen content, as expressed in Eq. (2).
Conversely, the growth respiration of vegetation is derived directly through empirical relationships. The model allocates 25% of the remaining organic matter for growth respiration after the subtraction of maintenance respiration from GPP. Therefore, the NPP is computed using Eq. (3), which integrates the calculated values of GPP, maintenance respiration, and growth respiration. This comprehensive approach ensures that the NPP estimations are reflective of the dynamic interplay between photosynthesis, respiration, and other metabolic processes within the ecosystem.
$${R}_{m}={R}_{{leaf}}+{R}_{{root}}+{R}_{{sapwood}}$$
(2)
$${\rm{NPP}}=0.75\times ({\rm{GPP}}-{R}_{m})$$
(3)
Where, Rm represent the carbon consumption from maintenance respiration, while \({R}_{{leaf}}\)、\({R}_{{root}}\) and \({R}_{{sapwood}}\) represent the carbon consumption from leaf, root, and sapwood respiration, respectively.
The original Fortran version of the LPJ-DGVM code is mainly divided into two parts: the core program encapsulated in LPJMAIN.F and the input-output interface program named LPJIO.F. The main program contains 22 specialized subroutines, while the input-output interface program contains 6 subroutines to manage data flow. The input data of the original code is monthly-resolved meteorological data, which is then interpolated to daily data using the daily function. To utilize daily-resolved meteorological data to drive the model, the LPJ-DGVM code has been refactored. The revised LPJ model can now produce parameters such as GPP, NPP, and NEP at daily resolution, thereby eliminating the need for interpolation and enhancing the temporal precision of the model’s inputs. At the same time, the spatial resolution of the model’s input data has been refined to a 5-km grid, and the control over soil and land use/cover data has been enhanced, improving the accuracy of the model’s simulation of surface cover. These enhancements are deliberate steps toward augmenting the model’s proficiency in conducting high-resolution simulations at the regional scale.
Vegetation Hydrological Balance Process Module and vegetation remote sensing data fusion
This specialized module within the LPJ-DGVM focuses on the water cycle processes specific to vegetation, including evaporation, transpiration, and interception of precipitation. It simulates how plants influence the local water balance by considering factors such as soil moisture content, plant water potential, and stomatal conductance. In the LPJ-DGVM model, the water cycle process simulation can be divided into a four-step process: (1) the computation of potential evapotranspiration, (2) the determination of actual evapotranspiration, (3) the estimation of soil moisture content, and (4) the quantification of runoff. From the standpoint of energy balance, the rate of evapotranspiration is intrinsically linked to temperature and net radiation. Leveraging this principle, the LPJ-DGVM employs Eq. (4) to estimate potential evapotranspiration.
$${E}_{pot}=\left(2.503\times {10}^{6}\times \frac{{e}^{\left(\frac{17.269T}{237.3+T}\right)}}{{(237.3+T)}^{2}\left)\right.}\right.\left/\left(\left(2.503\times {10}^{6}\times \frac{{e}^{\left(\frac{17.269T}{237.3+T}\right)}}{{(237.3+T)}^{2}}\right)+\gamma\right)\times {R}_{n}/{\rm{\lambda}}\right.$$
(4)
In the equation, the numerator represents the rate of change of saturated vapor pressure with temperature, \(\gamma\) represents the psychrometric constant, \({R}_{n}\) represents the net radiation, \({\rm{\lambda }}\) represents the latent heat of vaporization, and T is the temperature (°C).
Actual evapotranspiration is defined as the sum of vegetation transpiration and soil evaporation. The LPJ-DGVM considers various influential factors such as vegetation interception, soil moisture, and vegetation structure when calculating actual evapotranspiration. The calculation of vegetation interception uses the maximum water capacity of the vegetation canopy, \({S}_{I}\), which is a function of LAI and precipitation.
$${S}_{I}={{\rm{p}}}_{r}\times {\rm{LAI}}\times {\rm{i}}\times {{\rm{f}}}_{v}$$
(5)
Where, \({S}_{I}\) represents the maximum water capacity of the vegetation canopy and is defined as a function of LAI and precipitation. i is the interception coefficient, which is related to the type and structure of the biome. \({{\rm{f}}}_{v}\) represents the vegetation coverage. LAI is the LAI.
In this study, the LAI used is replaced by the remote-sensed LAI to replace the fixed value in the original model to better match the actual changes. This substitution from a fixed value to a dynamic, satellite-based LAI dataset allows for a more accurate reflection of the actual variations in vegetation density and structure across different landscapes and temporal scales.
The vegetation transpiration is taken to be the smaller of the environmental water demand and the vegetation soil water supply, as expressed in Eq. (6):
$${{\rm{E}}}_{T}={\rm{Min}}\left(\left({{\rm{E}}}_{\rm{max}}{{\rm{W}}}_{\!r}\right),\left((1-{\rm{w}}){{\rm{E}}}_{{pot}}\times \frac{{{\rm{\alpha }}}_{m}}{(1+{{\rm{g}}}_{m}/{{\rm{g}}}_{{pot}})}\right)\right)$$
(6)
Where, \({{\rm{W}}}_{r}\) is a coefficient calculated based on the soil moisture content in each layer and the distribution of vegetation roots. \({{\rm{E}}}_{\rm{max }}\) is a coefficient given based on specific vegetation. \({{\rm{\alpha }}}_{m}\) and \({{\rm{g}}}_{m}\) are two empirical coefficients that are given. \({{\rm{g}}}_{{pot}}\) is the optimal canopy conductance, which is calculated from the photosynthesis module. It is assumed that only the top 20 cm of soil has evaporation processes, and the calculation of soil evaporation is using Eq. (7), where \({{\rm{W}}}_{r20}\) is the relative soil moisture content of the top 20 cm layer. The total evapotranspiration is calculated.
$${{\rm{E}}}_{s}={{\rm{E}}}_{{pot}}\times \alpha \times {{\rm{W}}}_{r20}\left(1-{{\rm{f}}}_{v}\right)$$
(7)
In the LPJ-DGVM model, the soil is divided into two layers, and the changes in soil moisture content at different depths are calculated separately for the first and second layers through variables such as \({{\rm{E}}}_{T}\), \({{\rm{E}}}_{s}\), snowmelt, surface runoff, and groundwater runoff, as well as intercepted precipitation and soil infiltration.
The calculation of soil moisture content in the first and second layers is shown in Eq. (8), respectively.
$$\begin{array}{c}\Delta {w}_{1}=(p{r}_{t}+M-{\beta }_{1}\ast {E}_{T}-{E}_{S}-{P}_{1}-{R}_{1})/({w}_{\max }\ast {d}_{1})\\\Delta {w}_{2}=({P}_{1}-{\beta }_{2}\ast {E}_{T}-{P}_{2}-{R}_{2})/({w}_{\max }\ast {d}_{2})\end{array}$$
(8)
Where, \(p{r}_{t}\) is the rainfall value after subtracting vegetation interception, M is the amount of snowmelt, \({\beta }_{1}\) is the proportion of water absorbed from the first layer of soil during vegetation transpiration, \({\beta }_{2}\) is the proportion of water absorbed from the second layer of soil during vegetation transpiration, \({\beta }_{1}+{\beta }_{2}=1\). P1 and P2 is the infiltration amount of the first layer and the second layer, respectively. d1 and d2 are surface runoff and underground runoff.
If the calculated soil moisture content in the first layer exceeds the maximum soil moisture content of the first layer, the excess is considered as surface runoff. If the soil moisture content in the second layer exceeds the maximum soil moisture content of the second layer, the excess plus the soil water infiltration from the second layer is added as groundwater runoff.
The LPJ model has many intermediate variables, such as LAI, soil moisture, and soil temperature, which can be obtained directly through observations or calculated using other methods with higher precision. Data assimilation is an advanced technique that integrates new observation data into ongoing numerical simulations. This method considers the spatial and temporal distribution of data, as well as the potential errors in both the observational and background fields. By employing the data assimilation method to integrate direct observation data and remote sensing-derived data into the LPJ-DGVM, the accuracy of its simulation can be significantly improved. This paper uses MODIS-derived vegetation LAI to validate and optimize the LPJ model’s output LAI, iteratively refining the model’s algorithms. The data assimilation process in this paper employs the Ensemble Kalman Filter (EnKF) method, a widely recognized data assimilation algorithm particularly suited for nonlinear models. The EnKF method uses the Monte Carlo approach to solve the temporal process of probability density function during the simulation, which is composed of numerous model state ensembles. This probability density is composed of many model state collections, all of which are forward integrated over time from differential equations with random errors representing model uncertainties.
Other physiological and ecological process modules
Beyond the modules, the LPJ-DGVM also includes other physiological process modules. These include a module for simulating the allocation of photosynthetic products, a module for simulating the decomposition processes of litter and soil organic matter, a module for simulating vegetation death due to factors such as light competition and low growth efficiency, a module for simulating the generation of new individuals, and a module for simulating natural disturbances in ecosystems.
For herbaceous vegetation specifically, the allocation of photosynthetic products is simplified to occur primarily between leaves and roots. In the LPJ-DGVM, it is assumed that the decomposition rate of litter is 0.35 yr−1, the decomposition rate of rapidly decomposing soil organic matter is 0.03 yr−1, and the decomposition rate of slowly decomposing soil organic matter is 0.001 yr−1. Additionally, the LPJ-DGVM calculates the mortality of vegetation each simulated year and calculates the generation rate of new individuals each simulated year. The LPJ-DGVM also considers the impact of fire disturbances on ecosystems, which is generally believed to be primarily related to the amount of fuel and the moisture content of litter. Specific details will not be elaborated on here.
Parameterization of functional types of alpine vegetation (PFTs)
Due to the strong influence of the East Asian monsoon on our country, and the natural environment is relatively unique, the original model’s environmental limiting factors do not precisely capture the distribution of plant functional types in the Qinghai-Tibet Plateau. Based on in-depth analysis and synthesis of the vegetation distribution characteristics of the Qinghai-Tibet Plateau and the original 10 plant functional types of the LPJ model, and other studies33,34,35,36 of similar model parameterization schemes about Alpine Meadow Ecosystems, the new functional types have been developed and parameterized to better reflect the alpine vegetation distribution of the Qinghai-Tibet Plateau. New functional types, such as desert shrubs and cold-tolerant grass, have been introduced in adapting the model to the specific ecological conditions of the plateau (Table 2). Table 2 lists two added PFTs and their several physiological parameters in LPJ model. This customization has led to the establishment of an LPJ-DGVM model that is specifically tailored to simulate the alpine vegetation of the Qinghai-Tibet Plateau.