The main components of haze are SO₂, NO_X and particulate matter. The uptake of haze material by types of green space per unit area (Jin et al. 2005; Ye et al. 1998) (Table 1) and green space area were used to calculate the ecosystem service functions of haze absorption by green space (Han and Zhou 2015; Kuttler and Strassburger 1999).

The ecosystem service functions of haze absorption by green space include the absorption of SO₂, NO_X and respirable particulate matter. According to the various types of green space and the ecosystem service functions of haze absorption by each type of green space per unit area, the total ecosystem service functions of haze absorption by green space in China can be calculated from formula (1) (Han and Zhou 2015; Kuttler and Strassburger 1999).

$ESF=\sum _{i=1}^3\sum _{j=1}^3{A}_i{F}_{ij}$ (1)

where ESF is the total ecosystem service functions of haze absorption by green space; A_i is the area of green space type i; F_ij is the ecosystem service of absorbing haze component j by green space i per unit area; i is the green space type including forest cover, grass land and arable land; and j is the haze component including SO₂, NO_X and particulate matter.

Both the ecosystem itself and its spatial heterogeneity affect ecosystem service functions. Considering the ecological system, the quality of green space plays an important role in its function, and the vegetation coverage (normalized difference vegetation index (NDVI)) and net primary productivity (NPP) affect the corresponding service functions. The above ecosystem service functions calculation is only based on the land use area, without considering the impact of green space quality, so the results cannot reflect the true ecosystem service functions of haze absorption by green space. Using NDVI and NPP as evaluation indicators of green space quality and the correction coefficient to adjust the ecosystem service functions, the formula for the calculation is as follows (Gao et al. 2012):

$f_i=\frac{NDV{I}_i-NDV{I}_{\min }}{NDV{I}_{\max }-NDV{I}_{\min }}$ (2)

$Q_i=\left[\frac{NP{P}_i}{NP{P}_{\text{mean }}}+\frac{f_i}{f_{\text{mean }}}\right]/2$ (3)

$ES{F}^{\text{'}}=ESF\times Q_i$ (4)

where f_i and NPP_i are the NDVI and NPP of grid I, respectively; NPP_mean and f_mean are the mean NPP and NDVI values of various ecosystems in the study region, respectively; NDVI_max and NDVI_min are the maximum and minimum NDVI values for the entire growing season; Q_i is the green space quality coefficient; ESF is the ecosystem service functions before the green space quality correction; and ESF` is the ecosystem service functions after the green space quality correction.

To reflect the dependence of ecosystem service functions on the ecological functions index over time, the economic elasticity coefficient is selected to calculate the coefficient of sensitivity (formula (6)) (Kreuter et al. 2001).

$CS=\mid \frac{\left(ES{F}_j-ES{F}_i\right)/ES{F}_i}{\left(F_{jk}-F_{ik}\right)/F_{ik}}$ (5)

where ESF is the total ecosystem service functions; F is the functions coefficient; i and j are the initial and adjusted functions coefficients, respectively; k is the green space type; and CS is the coefficient of sensitivity. If CS > 1, the ESF for F is flexible, indicating that the total ecosystem service functions increase faster than the functions coefficient and that the proportion of the total ecosystem service functions and the functions coefficient are increasing. However, if CS < 1, the ESF for F is inelastic. CS = 1 represents complete elasticity; CS = 0 represents complete inelasticity. A higher ratio indicates that the elasticity of the ecosystem service functions index is more important.

Landscape pattern indices are used to describe the spatial organization of a landscape and provide a quantitative measure of the composition and spatial configuration of landscape structure. The interaction between landscape patterns and ecological processes as well as green space impacts haze absorption to different degrees. Based on previous research (Fang et al. 2014), we selected the landscape-level indices of patch density (PD), the interspersion and juxtaposition index (IJI), the area-weighted mean shape index (SHAPE_AM), and Shannon’s diversity index (SHDI) to study the relationship between landscape patterns and the ecosystem service functions of haze absorption by green space in China. Among these indices, SHAPE_AM was calculated by the formula from reference (Fang et al. 2014), and PD, IJI and SHDI were calculated by the following formulas.

$PD=\frac{1}{A}\sum _{i=1}^n{N}_i$ (6)

where PD is patch density; A is the total area of the landscape; N_i is the number of patches in landscape i; i is the landscape element; and n is the total number of patches in the landscape.

$IJI=\frac{-{\sum }_{i=1}^m{\sum }_{k=1+1}^m\left[\left(\frac{e_{ik}}{E}\right)·\ln \left(\frac{e_{ik}}{E}\right)\right]}{\ln (0.5[m-(m-1)\left]\right)}\times 100$ (7)

$IJI=\frac{-{\sum }_{i=1}^m{\sum }_{k=1+1}^m\left[\left(\frac{e_{ik}}{E}\right)·\ln \left(\frac{e_{ik}}{E}\right)\right]}{\ln (0.5[m-(m-1)\left]\right)}\times 100$ (8)

where m is the total number of landscape types; i and k are the numbers of patches of types i and k, respectively; e_ik is the total boundary length of the patch types between patch types i and k; E is the total boundary length of the landscape, including the background; and p_i is the perimeter of patch type i.

The ecosystem service functions of haze absorption by green space, including measures of the absorption of SO₂ and NO_X, dust retention and the total ecosystem service functions, were calculated for different provinces in China using a geographic information system (GIS). The calculations of landscape pattern indexes including PD, IJI, SHAPE_AM, and SHDI for provinces of China were performed in FRAGSTATS. Correlations between landscape patterns and the absorption of SO₂ and NO_X, dust retention and total ecosystem service functions were calculated as Pearson correlation coefficients as follows:

${\rho }_{x,y}=\frac{cov(X,Y)}{{\sigma }_X{\sigma }_Y}$ (9)

where cov (X, Y) represents the covariance between two variables, and σ_X and σ_Y refer to the variance of the two variables. The Pearson correlation coefficient is used to measure the correlation between two variables. The value of this coefficient falls between 1 and -1: 1 represents a full positive correlation of the variables; 0 indicates that the variables are independent; and -1 indicates a completely negative correlation.

A MODIS land cover classification product (mod12q1) was used for the land use data for China in 2001, 2004, 2007, 2010, 2013, 2016 and 2018. The spatial resolution of this product is 500 m, and land use is divided into arable land, forest cover, grass land, construction land, unused land and water bodies. Because the ecosystem service functions of haze absorption by water bodies are relatively small (Han and Zhou 2015), and few studies have been conducted on haze absorption by water bodies, it is difficult to obtain ecosystem service functions for the absorption of SO₂ and NO_X and dust retention by this land use type per unit area (Liu and Yu 2016), so the functions were not included as green space in this study. Therefore, green space in this study includes arable land, forest cover and grass land (Han and Zhou 2015; Kuttler and Strassburger 1999). Both NDVI and NPP are MODIS data products for China in 2001, 2004, 2007, 2010, 2013, 2016 and 2018 with a spatial resolution of 500 m, and in addition, the NPP data of NTSG (Numerical Terra-dynamic Simulation Group) was used as a supplement; the resolution of the data was 1 km × 1 km, and the annual NPP of the terrestrial ecosystem was obtained by using the NPP estimation model established by Biome-BGC and light energy utilization model. A dataset of the boundaries of the provinces in China was also included in this study.

As shown in Table 2, the total ecosystem service functions (benchmark values) of haze absorption by green space in China were 9000458.55 million Kg in 2001, 8784710.32 million Kg in 2004, 8900539.79 million Kg in 2007, 9179977.89 million Kg in 2010, 9145110.75 million Kg in 2013, 7761608.74 million Kg in 2016 and 7734526.75 million Kg in 2018. The ecosystem service functions of haze absorption by green space in China showed a trend of first increasing and then decreasing in 2001–2018, exhibiting an upward trend from 2001–2013 and increasing by 144652.20 million Kg (1.61%), primarily because the Chinese government invested 179 billion Yuan in a series of ecological restoration programs (Wang et al. 2007), including the Three North Shelterbelt Development Program, the Conversion from Cropland to Forest Program and the Natural Forest Protection Program, to restore degraded ecological environments and to foster stable and sustainable development. The area of forest cover increased by 6886085.77 ha (1.69%) between 2001 and 2013, but the ecosystem service functions fell by 1410584.00 million Kg from 2013–2018, a decrease of 15.42%, primarily because of adjustment of ecological land structure, the reduction of forest cover with high haze absorption ecological service function, the increase of grass land with low haze absorption function, and the reduction of arable land caused by the expansion of construction land. The ecosystem service functions decreased by 215748.23 million Kg from 2001–2004, a decrease of 2.40%, and the main reasons are a reduction in the area of forest cover and grass land. In contrast, the ecosystem service functions increased by 115,829.48 million Kg from 2004–2007, an increase of 1.32%, mainly due to the increase in forest cover and grass land and the decrease in arable land that is largely attributed to the Program for Conversion from Cropland to Forest and Grass Land in China. Additionally, the ecosystem service functions increased by 279438.10 million Kg from 2007–2010, an increase of 3.14%, primarily due to an increase in forest cover and a reduction in arable land and grass land. From 2010–2013, the ecosystem service functions decreased by 34,867.15 million Kg, a reduction of 0.38%, mainly due to the decrease in forest cover and grass land and the increase in arable land. In comparison, the ecosystem service functions fell by 1,383,502.01 million Kg from 2013–2016, a decrease of 15.13%, primarily because of the increase in grass land and decrease in arable land and forest cover. The ecosystem service functions decreased by 27081.99 million Kg from 2016–2018, a decrease of 0.35%, mainly attributable to a reduction in the area of forest cover and arable land.

The contributions to haze absorption by green spaces indicated that the types are very different (Table 2). The overall contribution of forest cover was the largest, and these proportions were 98.68%, 98.67%, 98.68%, 98.75%, 98.77%, 98.17% and 98.16% in 2001, 2004, 2007, 2010, 2013, 2016 and 2018, respectively. Grass land had the next largest contribution, accounting for 1.15%, 1.13%, 1.14%, 1.08%, 1.05%, 1.67% and 1.67% in 2001, 2004, 2007, 2010, 2013, 2016 and 2018, respectively. The total contribution of arable land was less than 1% from 2001–2018, mainly due to the large area of forest cover combined with the higher per-unit functions of respirable particulate matter and SO₂ in the haze, both of which resulted in a higher contribution to ecosystem service functions from the other types. In contrast, the relatively lower contributions from grass land and arable land were primarily due to the smaller per-unit functions of haze absorption.

The primary haze absorption ecological function by green space was primarily dust retention (Table 2), the functions of which accounted for 97.98%, 97.97%, 97.97%, 98.04%, 98.06%, 97.47% and 97.46% of the total functions in 2001, 2004, 2007, 2010, 2013, 2016 and 2018, respectively. Sulfur dioxide followed, accounting for 1.90%, 1.89%, 1.90%, 1.83%, 1.82%, 2.40% and 2.41% of the total functions in 2001, 2004, 2007, 2010, 2013, 2016 and 2018, respectively. The functions of NO_X absorption from 2001–2018 were less than 1% of the total, primarily due to the lower per-unit area function of the absorption of SO₂ and NO_X by various types of green space. However, the effect of dust retention was clear: the function of respirable particulate matter absorption by forest cover was especially high, indicating that green space plays an important role in dust removal and retention. Respirable particulate matter is the most important component of haze, so planning a reasonable amount of green space is conducive to reducing haze.

The ecosystem service functions (corrected value) of haze absorption by green space in China increased by 596086.46 million Kg (7.72%) from 2001–2013 (Table 3), while decreasing by 1810169.25 million Kg (21.76%) from 2013–2018: 7724215.34 million Kg in 2001, 7558663.89 million Kg in 2004, 7794138.04 million Kg in 2007, 8150709.06 million Kg in 2010, 8320301.79 million Kg in 2013, 6515690.35 million Kg in 2016 and 6510132.55 million Kg in 2018. In addition to the decrease of 165551.45 million Kg (2.14%) in 2001–2004, decrease of 1804611.45 million Kg (21.69%) in 2013–2016 and decrease of 5557.80 million Kg (0.09%) in 2016–2018, the functions from 2004–2007, 2007–2010 and 2010–2013 increased by 235474.15 million Kg (3.12%), 356571.02 million Kg (4.57%) and 169592.74 million Kg (2.08%), respectively.

The contribution rate of the various types of green space to haze absorption varied greatly. The contribution rate of forest cover was the largest, accounting for 98.80%, 98.86%, 98.89%, 98.97%, 99.03%, 98.29% and 98.22% of the total in 2001, 2004, 2007, 2010, 2013, 2016 and 2018, respectively. The second was grass land, which accounted for 0.96%, 0.86%, 0.84%, 0.79%, 0.72%, 1.58% and 1.65% of the total in 2001, 2004, 2007, 2010, 2013, 2016 and 2018, respectively. Arable land made the smallest contribution, accounting for 0.24%, 0.29%, 0.27%, 0.24%, 0.25%, 0.13% and 0.13% of the total in 2001, 2004, 2007, 2010, 2013, 2016 and 2018, respectively.

The haze absorption by green space was dominated by dust retention, the functions of which accounted for 98.09%, 98.15%, 98.18%, 98.26%, 98.32%, 98.30% and 98.29% of the total in 2001, 2004, 2007, 2010, 2013, 2016 and 2018, respectively. The function of SO₂ absorption accounted for 1.76%, 1.69%, 1.66%, 1.60%, 1.54%, 1.62% and 1.63% of the total, and the function of NO_X absorption accounted for 0.15%, 0.17%, 0.16%, 0.15%, 0.15%, 0.09% and 0.09% in 2001, 2004, 2007, 2010, 2013, 2016 and 2018, respectively.

As can be seen from Figure 1, the correction based on green space quality reduced the functions by 1276243.21 million Kg (14.18%), 1226046.43 million Kg (13.96%), 1106401.75 million Kg (12.43%), 1029268.84 million Kg (11.21%), 824808.96 million Kg (9.02%), 1245918.40 million Kg (16.05%) and 1224394.20 million Kg (15.83%) in 2001, 2004, 2007, 2010, 2013, 2016 and 2018, respectively, compared with the benchmark values of the ecosystem service of haze absorption by green space. However, in terms of general trends, the benchmark and corrected values of ecosystem services of haze absorption by green space in China increased by 144652.20 million Kg (1.61%) and 596086.46 million Kg (7.72%), respectively, from 2001–2013, while decreasing by 1410584.00 million Kg (15.42%) and 1810169.25 million Kg (21.76%), respectively, from 2013–2018, indicating that the ecosystem service functions based on green space quality differ greatly from the functions considering only the area green space. If only the green space area, and not the quality, is considered when evaluating the ecosystem service functions, the evaluation results will be too high. Nonetheless, the overall trends in the benchmark and corrected ecosystem service functions of haze absorption by green space in China are consistent during 2001–2013 and 2013–2018, showing first an increase and then a decreasing trend, indicating that the ecological restoration and conservation projects of the Chinese government have enhanced the ecosystem service functions of haze absorption by green space in 2001–2013, but the adjustment of ecological land structure and the expansion of construction land have led to a reduction in the ecosystem service functions of haze absorption in 2013–2018. The government should strengthen the restoration of forest vegetation with high haze absorption capacity and regulate the speed of urban expansion to improve the ability of haze absorption by ecological land.

Compared with the benchmark values, the contribution rates of the corrected value of haze absorption by forest cover increased by 0.11%, 0.19%, 0.21%, 0.21%, 0.26%, 0.12% and 0.06% in 2001, 2004, 2007, 2010, 2013, 2016 and 2018, respectively, while the contribution rates of the corrected value of haze absorption by grass land decreased by 0.18%, 0.27%, 0.30%, 0.29%, 0.33%, 0.09% and 0.02%, whereas those of arable land increased by 0.07%, 0.09%, 0.09%, 0.07% and 0.07% in 2001, 2004, 2007, 2010 and 2013 but decreased by 0.03% and 0.03% in 2016 and 2018, respectively.

The analysis of the ecosystem services functions of haze absorption by green space revealed that the corrected value of dust retention increased by 0.11%, 0.18%, 0.21%, 0.21%, 0.26%, 0.83% and 0.83% compared with the benchmark value in 2001, 2004, 2007, 2010, 2013, 2016 and 2018, respectively. The corrected value of SO₂ absorption decreased by 0.14%, 0.20%, 0.24%, 0.24%, 0.28%, 0.79% and 0.78% in 2001, 2004, 2007, 2010, 2013, 2016 and 2018, respectively. The corrected value of NO_X absorption in 2001, 2004, 2007, 2010 and 2013 increased by 0.03%, 0.03%, 0.03%, 0.03% and 0.02%, respectively, although the value decreased by 0.04% and 0.04% in 2016 and 2018, respectively. These results indicated that the benchmark and corrected values of the contribution rates of haze absorption by different types of green space and thus the ecosystem service functions are different, but all the functions exhibited a consistent trend. The contribution rates were ranked as forest cover, grass land and arable land, and the order of ecosystem service function was dust retention, SO₂ absorption, and NO_X absorption.

Figures 2–8 (benchmark values) show that the ecosystem service functions of haze absorption by green space had different spatial distributions in China from 2001–2018. In general, different ecosystem service functions had very different spatial distributions within the same year, while the spatial distribution of ecosystem service functions exhibited little difference between different years.

The maximum ecosystem service functions for the absorption of SO₂ (Fig. 2a), dust retention (Fig. 2c) and the total ecosystem services (Fig. 2d) for green space were 80459.35–100608.06 million Kg, 440884.21–8817697.00 million Kg and 69663.42–8882085.53 million Kg, respectively, in 2001. These services were primarily distributed in the northwestern, central-southern and northeastern regions, which is consistent with the spatial distributions of the different ecosystem service functions presented in Figures 3–8 (a, c, d) for 2004, 2007, 2010, 2013, 2016 and 2018, respectively. In contrast, the minimum ecosystem service functions for the absorption of NO_X (Fig. 2b) by green space were 0–964.80 million Kg in 2001, and high values for this service occurred mainly in the eastern and northeastern zones, which is inconsistent with the spatial distribution of NO_X absorption in Figures 3–8b for 2004, 2007, 2010, 2013, 2016 and 2018, respectively.

However, compared with the maximum and minimum values, intermediate ecosystem service functions for the absorption of SO₂ (Fig. 2a) by green space were 7719.35–27868.06 million Kg and 27868.06–80459.35 million Kg in 2001, and these functions were mainly distributed in the northwestern, southwestern, central, northeastern and eastern regions, which is consistent with the spatial distribution of the absorption of SO₂ in Figures 3–8a for 2004, 2007, 2010, 2013, 2016 and 2018, respectively. In addition, intermediate values for NO_X absorption (Fig. 2b) were 964.80–2266.90 million Kg and 2266.90–4024.24 million Kg in 2001, and these functions were mainly in the western, central-northern, central-southern, southern and southeastern regions, which is in accordance with the spatial distribution of the absorption of NO_X in Figures 3–8b for 2004, 2007, 2010, 2013, 2016 and 2018, respectively.

Intermediate ecosystem service functions for dust retention (Fig. 2c) by green space were 900.52–63259.00 million Kg and 63259.00–440884.21 million Kg in 2001, and these functions were mainly distributed in the northwestern, southwestern, central-northern and northeastern regions, which is consistent with the spatial distribution of dust retention in Figures 3–8 (c) for 2004, 2007, 2010, 2013, 2016 and 2018, respectively. Furthermore, intermediate values for total ecosystem services (Fig. 2d) were 8329.85–35420.81 million Kg and 35420.81–69663.42 million Kg in 2001, and these values were mainly in the northwestern, southwestern, southern, eastern and northeastern regions, which is in accordance with the spatial distribution of total ecosystem services in Figures 5–7d for 2010, 2013 and 2016, respectively, but is inconsistent with Figures 3–4 and 8d for 2004, 2007 and 2018, during which these functions were mainly distributed in the northwestern, southwestern, and central regions.

As shown in Figures 9–15 (corrected values), the maximum ecosystem service functions of SO₂ absorption (Fig. 9a), dust retention (Fig. 9c) and total ecosystem services (Fig. 9d) by green space were 53386.20–72412.13 million Kg, 5681893.88–7575858.50 million Kg and 5723384.25–7631179.00 million Kg in 2001, respectively, and were mainly in the southeastern, central-southern and southwestern areas, which is in accordance with the spatial distributions of the different ecosystem service functions presented in Figures 10–15 (a, c, d) for 2004, 2007, 2010, 2013, 2016 and 2018, respectively.

In contrast, the intermediate ecosystem service functions of SO₂ absorption (Fig. 9a) by green space were 0–10506.86 million Kg and 10506.86–53386.20 million Kg in 2001 and were mainly distributed in the southwestern, central-northern and northeastern regions, which is consistent with the spatial distribution of SO₂ absorption in Figures 10–15a for 2004, 2007, 2010, 2013, 2016 and 2018, respectively. Moreover, the intermediate ecosystem service functions of dust retention (Fig. 9c) by green space were 1893964.63–3787929.25 million Kg and 5681893.88–7575858.50 million Kg in 2001 and occurred mainly in the northeastern, central and southeastern areas, which is in accordance with the spatial distribution of dust retention in Figures 9–15c for 2004, 2007, 2010, 2013, 2016 and 2018, respectively. Furthermore, the minimum values of the total ecosystem service functions (Fig. 9d) by green space were 0–1907794.75 million Kg in 2001 and occurred mainly in the western, central-northern and southwestern zones, respectively, which is in accordance with the spatial distribution of the total ecosystem service functions (Figures 10–15d) in 2004, 2007, 2010, 2013, 2016 and 2018.

The results show that there was a great difference in the spatial distributions of the benchmark and corrected values of haze absorption by green space, and the spatial distributions of the maximum, intermediate and minimum ecosystem service function values were obviously different. However, the spatial distributions of the benchmark and corrected values also exhibited the same trend. In the same year, the spatial distribution of the ecosystem service functions of haze absorption by green space was very different, but in different years, the difference in the spatial distribution of the ecosystem service functions of haze absorption by green space exhibited little difference.

Figures 16–22 (benchmark values) show that the spatial distribution of ecosystem service functions and the proportion of haze absorption by green space differed in different provinces in China. Overall, different ecosystem service functions exhibited different spatial distributions in the same year or between different years. Some spatial distributions were quite different; others were more similar.

The ecosystem service functions of the absorption of SO₂ (Fig. 16a) and NO_X (Fig. 16b) by green space were 2.74–30586.00 million Kg and 0.16–1207.44 million Kg in 2001, respectively. The maximum and minimum values were primarily distributed in Xinjiang and in Shanghai, accounting for 17.88%, 0.02%, 10.98% and 0.14% of the total regional functions, respectively, The spatial distribution of these functions in 2001 was consistent with the absorption of SO₂ and NO_X in Fig. 17a and 17b but different from those shown in Figures 18–22a and 18–22b, which represent values for 2007, 2010, 2013, 2016 and 2018, respectively. Additionally, the functions of dust retention and total ecosystem services (Fig. 16c, d) for green space were 226.47–2729875.75 million Kg and 229.37–2761669.25 million Kg in 2001, respectively, and the maximum and minimum values for these services were primarily in Xinjiang and Shanghai, accounting for 30.97%, 0.01%, 30.70% and 0.01% of the regional totals, which is consistent with the spatial distribution of dust retention and total ecosystem services in Figures 17–22c and Figures 17–22d in 2004, 2007, 2010, 2013, 2016 and 2018, respectively. Most of the ecosystem service functions of haze absorption by green space were primarily from dust retention, which accounted for approximately 96% of the total. The functions for SO₂ absorption were the next highest, accounting for approximately 3% of the total, while NO_X accounted for approximately 1% (Fig. 16e) in 2001, which is consistent with the percentages of the ecosystem service functions of haze absorption by green space in Figures 17–22e in 2004, 2007, 2010, 2013, 2016 and 2018, respectively.

The ecosystem service functions of the absorption of SO₂ (Fig. 18a) by green space were 1.39–30872.46 million Kg in 2007, and the maximum and minimum values were mainly distributed in Xinjiang and Shanghai, accounting for 18.25% and 0.02% of the regional totals, respectively. This distribution is in accordance with the spatial distribution of the absorption of SO₂ in Figures 19–20a in 2010 and 2013, respectively. The ecosystem service functions of the absorption of NO_X (Fig. 18b) by green space was 0.10–1222.60 million Kg in 2007, and the maximum and minimum values were mainly in Xinjiang and Shanghai, accounting for 10.87% and 0.13% of the regional totals, respectively. This distribution is consistent with the spatial distribution of the absorption of NO_X in Figures 16b, 19b and Figures 21–22b in 2001, 2010, 2016 and 2018, respectively.

As shown in Figures 23–29 (corrected values), the ecosystem service functions of the absorption of SO₂ (Fig. 23a) by green space were 0.73–11,739.97 million Kg in 2001. The maximum and minimum values were primarily distributed in Yunnan and Shanghai, accounting for 8.62% and 0.03% of the total regional values, respectively. The spatial distribution of this value in 2001 was inconsistent with the absorption of SO₂ in Figures 24–29a, which represent values for 2004, 2007, 2010, 2013, 2016 and 2018, respectively. The ecosystem service functions of the absorption of NO_X (Fig. 23b) by green space were 0.04–1206.36 million Kg in 2001, and the maximum and minimum values were mainly distributed in Heilongjiang and Ningxia, accounting for 10.43% and 0.11% of the regional totals, respectively, which is consistent with the spatial distribution of absorption of NO_X in Figures 24–27b in 2004, 2007, 2010 and 2013, respectively, but inconsistent with Figures 28–29b in 2016 and 2018, respectively. The ecosystem service functions of the absorption of dust retention (Fig. 23c) and total ecosystem services (Fig. 23d) by green space were 58.62–930,837.56 million Kg and 59.39–943211.69 million Kg in 2001. The maximum and minimum values were primarily distributed in Yunnan and Shanghai, accounting for 12.29%, 0.01%, 12.21% and 0.01% of the total regional functions, respectively. The spatial distribution of this function in 2001 was consistent with the absorption of dust retention and total ecosystem services in Figures 24–27c and 24–27d, which represent functions for 2004, 2007, 2010 and 2013, respectively, but inconsistent with the Figures 28–29c and 28–29d in 2016 and 2018, respectively.

The results show that there was a great difference in the spatial distributions of the benchmark and corrected values of haze absorption by green space in different provinces in China, and the maximum and minimum of ecosystem service functions were obviously different. However, the spatial distributions of the benchmark and corrected values also exhibited the same trend. In the same year, the spatial distribution of the ecosystem service functions of haze absorption by green space was very different in different province, but in different years, the difference in the spatial distribution of the ecosystem service functions of haze absorption by green space exhibited little difference in different provinces.

The coefficient of sensitivity of the ecosystem service functions for different green space types was generally quite different from 2001–2018 (Table 4). The sensitivity coefficients for forest cover were elastic, while those of arable land and grass land were inelastic. The coefficients of sensitivity for forest cover were highest due to the large area of this cover type and the high ecosystem service functions coefficient for haze absorption by green space. The coefficients of sensitivity were 0.9868 in 2001, 2004 and 2007, 0.9875 in 2010, 0.9877 in 2013, 0.9817 in 2016 and 2018, respectively, and the change rates were ± 49.3424%, ± 49.3398%, ± 49.3405%, ± 49.3767%, ± 49.3832%, ± 49.0861% and ± 49.0842%, respectively. The coefficients of sensitivity for grass land were relatively small, with values of 0.0115 in 2001, 0.0113 in 2004, 0.0114 in 2007, 0.0108 in 2010, 0.0105 in 2013, 0.0167 in 2016 and 2018, and the change rates for grass land were ± 0.5733%, ± 0.5635%, ± 0.5702%, ± 0.5384%, ± 0.5268, ± 0.8365 and ± 0.8374%, respectively. The coefficients of sensitivity for arable land were the smallest due to the low ecosystem service functions coefficient of haze absorption by green space. The values of this coefficient were 0.0017 in 2001, 0.0019 in 2004, 0.0018 in 2007, 0.0017 in 2010, 0.0018 in 2013, 0.0015 in 2016, and 0.0016 in 2018, and the change rates were ± 0.0843%, ± 0.0967%, ± 0.0893%, ± 0.0850%, ± 0.0901%, ± 0.0774% and ± 0.0783%, respectively.

To quantitatively understand the relationship between land use patterns and ecosystem service functions, a correlation analysis was conducted (Table 5). There were significant correlations between many landscape pattern metrics and ecosystem service functions, which indicated that landscape patterns significantly affected ecosystem service functions. The correlation coefficients between PD and the ecosystem service functions of SO₂ absorption, NO_x absorption, dust retention and total ecosystem services exhibited significant negative relationships with correlation coefficients of -0.407, -0.511, -0.330 and -0.332, respectively. In contrast, the correlation coefficients between SHAPE_AM and the ecosystem service functions of SO₂ absorption, NO_x absorption, dust retention and total ecosystem services exhibited significant positive relationships with correlation coefficients of 0.650, 0.634, 0.568 and 0.570, respectively. These results indicate that PD and SHAPE_AM have important effects on different ecosystem service functions. In general, the larger the PD, the smaller the ecosystem service functions; the larger the value of SHAPE_AM, the greater the ecosystem service functions.

The correlation coefficients between IJI and the ecosystem service functions of SO₂ absorption, NO_x absorption, dust retention and total ecosystem services exhibited significant negative relationships with correlation coefficients of -0.606, -0.507, -0.449 and -0.452, respectively. The correlation coefficients between SHDI and the ecosystem service functions of SO₂ absorption, NO_x absorption, dust retention and total ecosystem services also exhibited significant negative relationships with correlation coefficients of -0.242, -0.316, -0.202 and -0.203, respectively. These results indicate that IJI and SHDI have important effects on different ecosystem service values. In general, the smaller the IJI and SHDI, the larger the ecosystem service functions.