Spatio-temporal assessing of natural vegetation regulation on SO2 absorption coupling ecosystem process model and OMI satellite data

Quantifying the contribution of natural ecosystems on air quality regulation can help to lay the foundation for ecological construction, and to promote the sustainable development of natural ecosystems. To identify the spatio-temporal dynamic changes of natural vegetation regulation on SO2 absorption and the underlying mechanism of these changes in Qinghai Province, an important ecological barrier and the unique natural ecosystems, the Biome-BGC model was improved to simulate the canopy conductance to SO2 and leaf area index (LAI) on the daily scale, and then the SO2 absorption by vegetation was estimated coupling SO2 concentration from satellite data. Our results showed that the annual average SO2 absorption of the natural ecosystems in Qinghai Province was 4.74 × 104 tons yr−1 from 2005 to 2018, accounting for about 40% of the total emissions. Spatially, the ecosystem service of SO2 absorption gradually decreased from southeast to northwest, and varied from 0 in Haixi state to 14.37 kg SO2 ha−1 yr−1 in Haibei state. The annual average SO2 absorption in unit area was 0.68 kg SO2 ha−1 yr−1, and significantly higher SO2 absorption was observed in summer with 0.45 kg SO2 ha−1 quarterly. The canopy conductance and LAI controlled by climate variables would be the dominant driving factors for the variation of SO2 absorption for natural ecosystems. The sensitivity analysis showed that SO2 concentration contributed more to the uncertainties of SO2 absorption than the conductance in this study. Our results could provide powerful supports for realistic eco-environmental policy and decision making.


Introduction
With the concern of environmental problems, especially the high increasing concentrations of air pollutants, the ecosystem services of air quality regulation are drawing more attention (Brander et al 2012, Baró et al 2014, Fisher et al 2014. Vegetations offer the ability to remove significant amounts of air pollutants, and consequently improve environmental quality and human health. Vegetations remove gaseous air pollutants primarily by absorption via leaf stomata, though some gases are removed by vegetation surface (Dochinger 1980, Lovett 1994, Nowak 1994. Sulfur dioxide (SO 2 ) is one of the most common air pollutants, mainly resulted from the burning of coal and oil, also emitted by some vehicles burning high-sulfur fuel, and by volcanic eruptions. It is severely harmful to the environment and human health, and the governance is difficult and high-cost. Quantifying SO 2 absorption provided by vegetation is vital for further realistic environmental policy and decision making.
SO 2 absorption measurements have been made using several methods, such as the micro meteorological method and chamber method for the short-term flux quantification, and the difference method based on sulfur content of plants for the long-term flux (Erisman andBaldocchi 1994, Ke et al 2021). Although direct measurement can obtain precise amounts of SO 2 absorption, it is expensive and often technically challenging. Currently, two common methods are used to assess regional scale SO 2 absorption by vegetation: empirical coefficient calculation and inferential modeling. The first method averages SO 2 absorption in unit area based on the difference measurements, and then multiplies by the total area of each biome to obtain the total SO 2 absorption (Costanza et al 1997, Nakamura et al 2005. Essentially, this is a static assessment method that ignores the vegetation dynamics in time and space, also the temporal-spatial heterogeneity of air pollutants (Nelson et al 2009). For inferential modeling, gaseous SO 2 absorption is inferred from the pollutants concentration in station, appropriate absorption rate and leaf area index (LAI), which is the amount of one-sided leaf area per ground area (Janhäll 2015, Manes et al 2016. Thus, rate-based modeling can reflect the variations of air quality regulation services caused by the heterogeneity of vegetation structure and growth, not just type change.
As to the inferential method, the determination of the canopy absorption rate and LAI is critical for estimation accuracy (Nowak et al 2006). Previous studies usually got these two variables separately, of which, the absorption rate was estimated based on the independent canopy conductance model or literature ranges, and the canopy LAI was obtained from ground measurement data or satellite product (e.g. MODIS LAI). However, they had been almost no consideration for the mutual coupling and feedback of these two variables, which may bring in some estimation errors. For example, soil moisture is the critical intermediate factor coupling vegetation and environment variables, and the essential limit for canopy conductance (Roberts et al 1984, Harris et al 2004, Kumagai et al 2004, but often ignored. As an important ecological barrier and the unique natural ecosystems, the Qinghai Province provides vital ecosystem services from regional to global scales, such as carbon sequestration, water supply, biodiversity conservation, air quality regulation, flood regulation, cultural services, etc (Dong and Sherman 2015, Xu et al 2019, Liu et al 2020. Particularly, the high-cold ecosystems in Qinghai Plateau play an important role in climate regulation and air quality regulation in Asia. However, previous studies for air quality regulation services mainly focused on urban ecosystem (Nowak et al 2006, Selmi et al 2016, Parsa et al 2019, little research for natural ecosystems specifically. In terms of area proportion globally, the green area of natural ecosystems is far more than that of urban green space. Besides, due to the dynamic characteristics of the atmosphere, substantial amounts of gaseous pollutants are transported to the areas distant from the source (Abou Rafee et al 2017). Natural ecosystems also play an important role in reducing pollutant concentrations at regional, continental, and global scales. Therefore, more efforts are needed to quantify the contribution of natural ecosystems to air purification, which can help to lay the foundation for ecological construction, and to promote the sustainable development of natural ecosystems (Naidoo et al 2008, Barnaud andAntona 2014).
In this study, we took Qinghai Province for instance, and combined ecosystem process model with satellite-derived air pollutant data for quantifying natural vegetation regulation on SO 2 absorption. An improved Biome-BGC model that simulates the canopy LAI and daily absorption rate synchronously, also ozone monitoring instrument (OMI) satellite air pollutant data were adopted to explore the spatiotemporal characteristics of air quality regulation. The main purposes of this research were: (a) to quantify the SO 2 absorption service of natural ecosystems in Qinghai Province from 2005 to 2018; (b) to identify the spatial dynamic changes of services' quantity, and the absorption capacity among different biomes; (c) to explore the temporal variability of SO 2 absorption service, and the underlying influence mechanism.

Study area
The study was implemented in Qinghai Province, an important part of the Qinghai-Tibet Plateau ecological barrier. It is located in the south-central part of northwestern China (figure 1), and the northeastern of the Qinghai-Tibet Plateau. The Qinghai Province covers an area of 71.75 × 10 4 km 2 , extending from 89 • 35 ′ to 103 • 04 ′ E in longitude and from 31 • 39 ′ N to 39 • 19 ′ N in latitude. The predominant climate in this region is plateau continental climate, with a cool, short summer and a cold, long and dry winter. The elevation ranges from 1650 to 6860 m, the mean annual temperature is 0.88 • C, and the annual rainfall is about 350 mm. The vegetation type is diverse, and the alpine grassland is the main vegetation type in the study area (figure 1). From southeast to northwest, the production and species gradually reduce, in turn presenting forest, grassland and desert ecosystems. Wetland is less distributed, and the distribution of farmland is mainly concentrated in the low-altitude areas of the east.

Biome-BGC model
In this study, an improved Biome-BGC model was developed to quantify the service of SO 2 absorption in Qinghai Province. Biome-BGC model is a mechanistic biogeochemical model designed to simulate the storage and flux of water, carbon, and nitrogen within terrestrial ecosystems (Running and Hunt 1993  This model uses meteorological data (maximum and minimum temperatures, precipitation, solar radiation and vapor pressure deficit, etc), site-specific data (effective soil depth, soil texture, site elevation, latitude, nitrogen deposition and biological nitrogen fixation), and ecophysiological data (e.g. maximum stomatal conductance, fraction of leaf N in Rubisco, canopy specific leaf area, etc) to simulate the biogeochemical processes of the given biome at daily timestep (Hidy et al 2012).
In the Biome-BGC model, the carbon flux block includes photosynthesis, respiration, allocation, decomposition and mortality processes. Photosynthesis is a process in which CO 2 , water molecules and solar energy are combined to generate simple sugars by vegetation. This reaction is calculated through Farquhar's photosynthesis routine (Farquhar et al 1980), characterized using four equations. First, the photosynthetic rate follows the law of diffusion (equation (1)), meanwhile, it is assumed to be constrained by the carboxylation rate of photosynthesis reaction and the electron transport of RuBP regeneration (equation (2)). (1) where, A is the net photosynthetic rate, C a is the atmospheric CO 2 concentration, C i is the intercellular CO 2 concentration, g CO2 is the canopy conductance to CO 2 . A c is the RuBP carboxylase (Rubisco) limited rate of carboxylation, A j is the rate of carboxylation allowed by the capacity to regenerate RuBP, R d is the dark respiration (daytime leaf maintenance respiration); V cmax is the maximum rate of carboxylation, Γ * is the CO 2 compensation point, K c and K o are the kinetic constants for rubisco carboxylation and oxygenation; J max is the maximum rate of electron transport. The total CO 2 assimilation can be obtained based on the net photosynthetic rate and canopy LAI. The conductance in unit leaf area is calculated by stomatal, boundary layer, cuticle level conductance and mesophyll conductance (Baldocchi et al 1987). Since numerous field studies have shown that mesophyll resistance for CO 2 and SO 2 are small (Nowak et al 2006, Wesely 2007, thus we assumed it to be zero for CO 2 and SO 2 . The specific equation is as follows: where g bl , g s and g c represent the boundary layer conductance, stomatal conductance and cuticle level conductance, respectively. Stomatal conductance was calculated using the modified Jarvis model (Jarvis 1976, Damour et al 2010. From equation (6), stomatal conductance (g s ) is the product of the response functions to environment factors (soil water status, light, air humidity, temperature and atmospheric CO 2 concentration), each function being generally determined by boundary line analysis (Chambers et al 1985, Yu et al 2017. The hypothesis behind this approach is that the response to each environmental factor is independent from the others. where, g max represents the maximum stomatal conductance, f (T), f (VPD), f (C a ), f (PPFD) and f (Ψ ) are the response functions of air temperature (T), vapor pressure deficit (VPD), atmospheric CO 2 concentration (C a ), photosynthetic photon flux density (PPFD) and soil water potential (Ψ ), respectively. The Biome-BGC model used in this study, considers not only the influence of meteorological factors and atmospheric CO 2 concentration on canopy conductance, but also the responses and feedbacks of the terrestrial ecosystems process to canopy conductance, mainly reflected in the variation of soil moisture. Hence, more realistic and accurate absorption rate and LAI may be obtained, using a complete ecosystem process model.

SO 2 absorption module
As the absorption of SO 2 is directly related to the CO 2 uptake by photosynthesis and the H 2 O release by transpiration, according to Graham's law of diffusion (Slovik et al 1996), the canopy conductance to SO 2 in unit leaf area is proportional to CO 2 (and H 2 O): where g SO 2 is the canopy conductance to SO 2 , g CO2 is the canopy conductance to CO 2 , M CO2 and M SO2 refer to the molecular mass of CO 2 and SO 2 , respectively. Here, a simplified SO 2 absorption module was specifically developed based on the original Biome-BGC model. Similar to CO 2 , SO 2 enters into the leaves primarily through the stomata following the same diffusion pathway as CO 2 (Abrol and Ahmad 2003). Therefore, the amount of SO 2 absorbed by vegetation is determined by the concentration in the air, the concentration in the pore cavity, the conductance of the gas into the leaves, as well as the LAI (Hällgren et al 1982, Slovik et al 1996, Hu et al 2016: (8) where F SO2 is the total absorption of SO 2 by vegetation, A SO2 is the net absorption rate per unit leaf area, S a is the atmospheric SO 2 concentration, S i is the intercellular SO 2 concentration. The intercellular SO 2 concentration is assumed to be in equilibrium with the SO 2 concentration in the aqueous phase within the leaf, and this equilibrium is determined by the solubility of SO 2 in water (Kropff 1989). SO 2 entering leaves will dissolve rapidly in water and form into HSO 3 -, SO 3 2and H + , and then these highly toxic compounds are oxidized to sulfate, which is less toxic than sulfide. Sulfate can also be incorporated into organic compounds by the sulfate assimilation process. When the atmospheric SO 2 concentration is no more than the saturation value, SO 2 entering leaves is dissolved and assimilated rapidly, thus the intercellular concentration of SO 2 is approximately zero.
With the increase of atmospheric SO 2 , the SO 2 absorption may exceed the assimilation requirement of vegetations. In this case, it is not comprehensive for just considering the absorption. In this study, the SO 2 assimilation was calculated based on the CO 2 assimilation, as well as the ratio of carbon and sulfur in plant (Tabatabai 1984). When SO 2 absorption is greater than the maximum assimilation, we adopted the assimilation as the SO 2 absorption service provided by vegetation. As the atmospheric SO 2 concentration continues to rise, excessive sulfide gradually accumulates in the plant, which may reduce stomatal conductance and even induce stomatal closure (Omasa et al 1985, Kropff 1989. Nevertheless, the atmospheric SO 2 concentration in Qinghai Province has not reached the threshold (Biscoe et al 1973) that stomata respond, thus stomatal behavior affected by SO 2 was not considered here.
For the output control of the Biome-BGC model, we added the variables of the sulfur fluxes at daily, monthly, and annual time scales. By the above methods, the actual daily absorption of SO 2 was obtained, and then the monthly and annual absorption of SO 2 can be estimated through accumulation.

Meteorological dataset
The daily meteorological data for model-driven includes the daily maximum temperature, daily minimum temperature, daily average temperature, daily precipitation, daily average VPD, average solar shortwave radiation, and day length. They were obtained from ERA5-Land reanalysis dataset (www.ecmwf.int/en/forecasts/datasets/reanalysisdatasets/era5). ERA5-Land is a reanalysis dataset providing a consistent view of the evolution of land variables over several decades with an enhanced resolution of 0.1 • × 0.1 • . In order to simulate the diversity of ecosystems more accurately, the resolution of this dataset was downscaled to 1 km by the random forest algorithm (Holden et al 2011). Except VPD and day length, other variables are available in the downscaled dataset. VPD was obtained based on the relative humidity and temperature, and day length was calculated using the number of days and latitude. In addition, the planetary boundary layer height (PBLH) is extracted from ERA5-Land reanalysis, to calculate the surface SO 2 concentration.

Land use dataset
Based on the multi-year national resource continuous inventory data of China, including forest resource inventory, wetland resource survey, and land resources survey, we obtained the base map of the ecosystem types. Considering the overlapping regions of different data sources, auxiliary data and expert knowledge were used to redefine the type of ecosystem. The multi-source data as an auxiliary consists of the land use classification data from remote sensing, a 1:1000 000 scale vegetation distribution map of China, and sample plot surveys. The map of ecosystem types from multi-data fusion is presented in figure 1.

Air pollutant dataset
The OMI is an instrument boarded on the Aura satellite of the NASA Earth Observation System, that monitors atmospheric O 3 , NO 2 , SO 2 , and other trace gases (Levelt et al 2006). The OMI provides SO 2 column concentrations at four different heights: the PBL, the bottom of the troposphere, the middle of the troposphere, and the top of the troposphere to the stratosphere (Krotkov et al 2008, Wang et al 2015. In this study, we chose the daily Level-3 OMI PBL SO 2 product (OMSO2 v1.2.0) from 1 January 2005 to 31 December 2018, with the spatial resolution of 0.25 • × 0.25 • , which is available on NASA website (http://disc.sci.gsfc.nasa.gov/Aura/data-hold ings/OMI/omso2_v003.shtml).
The daily Level-3 data is estimated from the OMI SO 2 L2 product using the principal component analysis algorithm, and contains SO 2 vertical column densities (VCDs) in Dobson units (1 DU = 2.69 × 10 16 molecules cm −2 ). The quality of daily SO 2 satellite data was generally poor, and there were a large number of missing values and noisy points. Therefore, preprocessing and quality control were required to eliminate the poor-quality pixels in this study. The screening criteria were as follows: (a) radiative cloud fraction less than 0.2, and (b) solar zenith angle less than 70 • . Furthermore, spatial interpolation and spatial filtering were conducted to improve the range and quality of the available data. The data was processed to monthly mean values at a spatial resolution of 0.1 • × 0.1 • .
The PBL SO 2 data presents the column concentration of PBL, but the ground level concentration is required for estimating SO 2 absorption by vegetation. Previous studies have demonstrated that the satellite is an effective way to monitor the monthly-averaged characteristics of the surface SO 2 (Knepp et al 2015, Zheng et al 2018. Hence, this study employed the empirical method to convert the column concentration to the ground level concentration. Specifically, considering SO 2 varies weakly with heights in the PBL (Zhang et al 2019), the correlation between PBL scaled column SO 2 (VCD SO 2 /PBL) and groundbased measurements of SO 2 concentration without raining was built, and then the gridded surface SO 2 concentration could be obtained.

Vegetation contribution to SO 2 absorption
From 2005 to 2018, the annual average SO 2 absorption in Qinghai Province was 4.74 × 10 4 tons. Due to the largest area proportion, the contribution of grassland for SO 2 absorption was the largest in Qinghai Province. The annual average SO 2 absorption of grassland ecosystems was 2.67 × 10 4 tons (figure 2(a)), accounting for 56.18% of the total, and shrubland ecosystems made up 14.88% of the total absorption amount of SO 2 . Although the desert ecosystems made up 20.77% of the total areas, contributed only 2.42% to the total SO 2 absorption in Qinghai Province.
According to the data from the China Environment Statistical Yearbook, over the last 10 years, 13.90 × 10 4 tons of SO 2 was emitted annually for Qinghai Province (National Bureau of Statistics (NBS) and Ministry of Environmental Protection (MEP) 2005-2018). Here, we estimated the annual average SO 2 removal by vegetation absorption in Qinghai Province was 4.74 × 10 4 tons, accounting for about 40% of the total emissions. The ecosystem services of air purification refer to the benefits to humans due to the ecosystem self-purification ability, specifically, resulting from the reduction in the cost of manual decontamination of air pollutants. Without this part of the ecosystem purification function, we would have to invest more millions of dollars in cleaning up these pollutants (Nowak et al 2006).

Spatial distribution of SO 2 absorption
In terms of spatial distribution, the ecosystem service of SO 2 absorption gradually decreased from southeast to northwest, and varied from 0 in Haixi state to 14.37 kg SO 2 ha −1 yr −1 in Haibei state (figure 3). The regions with the SO 2 absorption greater than 2.50 kg SO 2 ha −1 yr −1 are mainly distributed in the eastern part of Qinghai Province. The annual average SO 2 absorption within Huangnan, Xining, and Haidong city were 2.09, 1.91, and 1.76 kg SO 2 ha −1 yr −1 , respectively. As the landscape has virtually no vegetation in the western desert region, the SO 2 absorption by vegetation was nearly zero, mainly distributed in Haixi state.
The SO 2 absorption is controlled by the SO 2 concentration and land cover patterns in the spatial distribution. The majority of population is concentrated in the eastern regions of Qinghai Province, mainly in Xining and Haidong districts, which was highly accordant with the spatial distribution of SO 2 concentration ( figure 4(a)). Additionally, from southeast to northwest, the production and plant species gradually reduce, in turn presenting forest, grassland and desert ecosystems (figure 1). As the most important determinant for total SO 2 absorption in Qinghai Province, SO 2 absorption capacity varied quite a lot among different vegetation types.
The annual average SO 2 absorption of Qinghai Province in unit area was 0.71 kg SO 2 ha −1 yr −1 over the past decade, and the highest SO 2 absorption service in unit area occurred in forest ecosystems (2.08 kg SO 2 ha −1 yr −1 ), slightly above shrubland  ecosystems (1.51 kg SO 2 ha −1 yr −1 ). The annual average SO 2 absorption service in unit area of grassland ecosystems was 0.85 kg SO 2 ha −1 yr −1 , around the same as wetland ecosystems. The greater capacity of SO 2 absorption in forest ecosystems mainly attributes to the high LAI. The SO 2 absorption capacity of vegetation was consistent with the photosynthetic capacity for the same diffusion pathway as CO 2 .

Interannual and seasonal variation of SO 2 absorption
In the past 10 years, the annual SO 2 absorption of Qinghai Province displayed a slight increase with fluctuation ( figure 2(b)). The SO 2 absorption service was the lowest in 2012, with the absorption in unit area of 0.62 kg SO 2 ha −1 yr −1 , and the highest was 0.78 kg SO 2 ha −1 yr −1 in 2011. All natural ecosystems were showing the rises of SO 2 absorption, but the time of the peak absorption varied among different ecosystems. For example, the SO 2 absorption services peaked in 2007 and 2018 for forest and grassland ecosystems, respectively. Moreover, the absorption valley of wetland ecosystems occurred in 2011, different from other ecosystems, but same with the total. Figure 5 shows the monthly mean temperature, precipitation, SO 2 concentration and SO 2 absorption from 2005 to 2018. For Qinghai Province, it presented a higher SO 2 concentration in winter and spring ( figure 5(b)), and relatively lower concentration in summer and autumn, which was mainly determined by two aspects. The low vegetation cover in winter greatly reduced the absorption of air pollutants. Meanwhile, according to the Statistical Review of China energy, the coal consumptions of Qinghai  Province in the 4th and 1st quarters are obviously higher than that in others, thus the combustion of sulfur-containing fossil fuels in winter further accelerates the SO 2 emissions to the air. Things were the opposite for the summer, the high vegetation cover in the growing season played a positive role in reducing SO 2 concentration.
Among the four seasons, the SO 2 absorption in summer (June-August) was observed significantly higher than that in others, with a total amount of 0.45 kg SO 2 ha −1 . The SO 2 absorption reached the maximum in August, while SO 2 concentration showed a reverse change trend. In winter (December-February), the SO 2 absorption service provided by vegetation was close to 0. Moreover, the SO 2 absorption services were 0.07 and 0.21 kg SO 2 ha −1 in spring (March-May) and autumn (September-November), respectively. For forest and shrubland ecosystems, the SO 2 absorption services in spring were much greater than that in others. The dominant species in Qinghai Province are cypress and spruce, belonging to temperate evergreen coniferous forest, thus forest and shrubland ecosystems own longer growing seasons. The SO 2 absorption service of desert ecosystems was less than 0.2 kg SO 2 ha −1 in the four seasons, and near to zero except in summer. The SO 2 absorption service of cropland ecosystems in summer was 1.60 kg SO 2 ha −1 , the highest among all ecosystems because the growing season is more concentrated.

Comparison with previous research
The rate-based method takes the concentration of pollutants into account, and usually calculates air pollutants absorption combining pollutants concentration, absorption rate and canopy LAI, which can reflect the dynamic service of SO 2 absorption. Hence, this method was more commonly applied to assess the air-quality regulation service, especially in urban green land and forest ecosystems. For example, Nowak et al (2006) estimated average leaf-on daytime dry deposition velocities ranged from 0.38 to 0.69 cm s −1 for SO 2 , and absorption service of urban trees ranged from 1 to 35 kg SO 2 ha −1 yr −1 among different cities in the US. Baró et al (2014) analyzed the seasonal differences of SO 2 absorption for the urban forests in Barcelona, showing that spring and summer were the seasons with higher absorption rate in average (the percent of uptake during the two seasons was 70.46% for SO 2 ). Table 1 summarizes some research for estimating SO 2 absorption by the ratebased method.
In this study, based on the improved Biome-BGC model, the SO 2 absorption by vegetation was 0.71 kg SO 2 ha −1 yr −1 in Qinghai Province, 2.08 and 1.51 kg SO 2 ha −1 yr −1 for forest and shrubland ecosystems respectively, which was relatively consistent with the research results in Barcelona and Tabriz, but lower than the absorption capacity of urban green land in US and Korean. Considering high pollutant concentration and vegetation cover of urban green land in US and Korean, the SO 2 absorption service provided by vegetation in the urban areas was greater than that of natural ecosystems, thus our results were deemed reasonable.
By contrast, the empirical coefficient method ignores the spatio-temporal heterogeneity of vegetation canopy and air pollutants, thus just presenting the total amount of the SO 2 absorption, lacking details in different time and space. In addition, the unit coefficient from difference measurements may be overvalued, because the sulfur uptake from soil is in. Therefore, this method usually overestimates gaseous SO 2 absorption by vegetation. For example, according to the unit coefficient from China's Biodiversity Situation Research Report (Zhao et al 2004), the SO 2 absorptions by vegetation were estimated to reach 0.73, 2.97 and 14.5 × 10 4 tons for forest, shrubland and grass ecosystems in Qinghai Province, apparently higher than the results in this study (figure 2).

Underlying factors for SO 2 absorption variations
The SO 2 concentration, canopy conductance in unit leaf area and canopy LAI lead to the fluctuation of SO 2 absorption together (equation (8)). From figure 4(b), SO 2 concentration in Qinghai Province was ascending as time passed in the last 10 years, while no obvious change of SO 2 absorption was found. Seasonally, although the high SO 2 concentration in winter and spring, the low vegetation cover still limits the absorption of SO 2 , thus the absorption was concentrated mostly in June-September (figure 5). To eliminate the impact of canopy LAI, we compared the SO 2 concentration and absorption rate per unit leaf area A SO2 in June-September. The SO 2 concentration in June is much higher than that in August, but A SO2 in August is significantly higher than that in June (figure 6(c)). This suggested no significant correlation between the SO 2 absorption and concentration in Qinghai province, even with the SO 2 emissions increased by human activity.
Here, we selected the SO 2 concentration, LAI, and SO 2 absorption from June to September during 2005-2018 to further analyze the dominant factor for SO 2 absorption, as shown in figure 6. Correlation analysis showed that canopy LAI by climate variables would be the dominant driver for the variation of SO 2 absorption for natural ecosystems, where pollutant concentration is relatively low. Likewise, the canopy conductance g so2 is also driven by climatic factors, and analysis is conducted in the next paragraph. Figure 6 also presents an extremely high absorption in June, 2011. This is caused by the accidental high concentration of SO 2 resulted from the volcano eruption of Nabro (Clarisse et al 2014), which brings in a short-term increase of SO 2 concentration regionally. In conclusion, in a normal situation when the pollutant concentration is within a certain threshold, climate variables are the most important determinants of vegetation SO 2 absorption for natural ecosystems.
For the alpine ecosystem, temperature and precipitation are two key factors for vegetation growth, and control carbon and sulfur uptake of ecosystem directly or indirectly through affecting plant stomatal activities (Ma et al 2019). The seasonal variation of SO 2 absorption service agreed with temperature and precipitation. In spring, the temperature was the dominant factor controlling SO 2 absorption. It might be due to that increasing spring temperature could stimulate photosynthetic enzyme activities and improve the availability of soil water (Bergh andLinder 1999, Welp et al 2007), and then ignite CO 2 and SO 2 absorption indirectly through its impacts on stomata. In summer, rising precipitation enhances soil water content and promotes vegetation transpiration further, which would lead to stomata respond and strengthen SO 2 absorption indirectly (Saliendra et al 1995).

Uncertainties and limitations
The uncertainties mainly result from the estimates of the conductance and satellite-derived surface concentration. Considering the complexity of the Biome-BGC model, the quantification of the pollutants absorption rate remains a degree of uncertainty. In addition, the errors in the OMI observations and empirical model can lead to uncertainties in surface concentrations of SO 2 . Here, we made a simple sensitivity analysis to quantify the uncertainties of these two factors roughly. Specifically, the estimation errors were recorded with the percentage increase and decrease of inputs (±5%, ±10%, ±15%, ±20%). The results showed that the SO 2 concentration reached higher estimation errors than the conductance, with the same percentage changes. Moreover, with the increasing error on SO 2 concentration, the uncertainties of absorption estimate increased more. Therefore, accurate SO 2 concentration near the surface is the key to improving the estimate accuracy of SO 2 absorption by vegetation on a regional scale.
At present, the sulfur cycle and its effects on stomatal conductance are not considered in the Biome-BGC model. When atmospheric SO 2 concentration is high, SO 2 absorbed by vegetation could not be completely assimilated, and even out from vegetation. In order to minimize the influence of the errors from model structure, a modified algorithm was also proposed. We adopted the assimilation as SO 2 absorption service provided by vegetation when the SO 2 absorption is greater than the assimilation, which was approximately obtained based on the stoichiometry method. In future research, the specific model for the fluxes of SO 2 into leaves and assimilation is required for assessing air purification service accurately.

Conclusions
In this study, we carried out the assessment of natural vegetation regulation service on SO 2 absorption in Qinghai Province, based on biogeochemical modeling and satellite-derived data. This was achieved with an improved Biome-BGC model, which is a targeted tool to simulate the storage and flux of water, carbon, and nitrogen within terrestrial ecosystems. The coupling of satellite air pollutant data with the modeling of ecosystem processes is a powerful way to identify the spatio-temporal dynamics of air purification services and explore the underlying mechanism, that can also be used to assess other services provided by natural ecosystems, such as climate change mitigation, water regulation, etc, thereby fully exploiting the multi-functional role of the ecosystem.
Our research showed that the annual average SO 2 absorption service of the natural ecosystems in Qinghai Province was 4.74 × 10 4 tons, accounting for about 40% of the total emissions. The contribution of natural ecosystems to air purification is significant enough to maintain ecosystem balance and reduce the cost of air pollution treatment. Air pollutant absorption by natural vegetation is one of the key components of natural assets, and its accurate assessment can also provide strong supports for ecological construction and environmental management.

Data availability statement
No new data were created or analyzed in this study.