Global relative ecosystem service budget mapping using the Google Earth Engine and land cover datasets

Ecosystem service mapping (ESM) studies are receiving increasing attention due to the imbalance between the supply of and demand for ecosystem services (ES). Global scale ESM is still scarce, but the high computing power of the Google Earth Engine (GEE) cloud platform significantly increases the efficiency. Based on global-scale land cover datasets and the GEE, an ES matrix model based-expert is constructed in this paper to map the ES supply, demand, and relative budgets. The net primary productivity (NPP), enhanced vegetation index (EVI), nighttime light (NTL), and world population (Pop) were acquired, and the NPP and EVI and the NTL and Pop datasets were used to revise the supply of and demand for ESs, respectively. We discovered that the ES supply capacity exhibits a double-peaked distribution with latitude, and the peaks are located at the equator and 50° N. The global ESs have a high spatial heterogeneity and the global supply of ESs is 2.405 times higher than the demand; however, the demand exhibits an increasing trend of about 3.36% per decade, and only southern Asia has more ES demand than supply. The imbalance between the ES supply and demand produced a push-pull effect, that is, it forced humans to move closer to the ES surplus regions (ESSRs) and farther away from the ES deficit regions (ESDRs), and the destruction of the ecological environment promoted this phenomenon. The global terrestrial area is divided into eight ES sub-regions, and targeted land management, urban planning, and environmental remediation policies are proposed.


Introduction
Ecosystem services (ESs) are defined as the benefits that ecosystems directly or indirectly provide to humans (Costanza et al 1997a, b, Ouyang andWang, 1999). ES mapping is the process of quantitatively describing the spatial distribution of ESs at spatiotemporal scales under the combined influence of natural and anthropogenic factors (Verburg 2014, Essen 2015. The ES mapping results are a crucial component in the formulation and implementation of ecological and environmental protection planning and management decisions. The ES flow caused by the imbalance between the supply of and demand for ESs provides the specific and quantifiable spatiotemporal characteristics of biophysical information for decision-makers to construct optimal policies. Since research on ecosystem service mapping (ESM) promotes sustainable development and the benefits derived from natural ecosystems, it has received increasing attention (Fu et al 2013, Mononen et al 2014, Drius et al 2018. With the rapid research of ES, especially ESM, the number of studies conducted at diverse scales (e.g., the landscape, regional, and national scales) (Nahuelhual et al 2013, Paracchini et al 2014, Haase et al 2016 is growing rapidly, however, relatively few studies have been conducted on the global scale (Verburg et al 2013). The fact that the Earth is a macro ecosystem with carbon sequestration, climate regulation, air circulation, and pollination cannot be ignored and needs to be studied at the global scale (Arrigo et al 2020, Nankaya et al 2021. Thus, it is urgent and necessary that global scale ESM be conducted to determine the ES supply, demand, and relative budgets. ES supply and demand are the primary components of ESM (Syrbe and Grunewald, 2017), and they can be used visualize the spatial imbalance between the materials provided by ecosystems and the human demand for these materials (Xiao et al 2016, Guan et al 2020. The ES supply and demand correspond to the ecoenvironment and socio-economic components (Yang et al 2012), respectively. Mapping the ES supply and demand is a crucial and straightforward approach to the harmonization of both components. The destruction of natural ecosystems and the dramatic reduction of the ES supply capacity are primarily due to human activities (e.g., deforestation, urban sprawl, and carbon emissions) (Yang et al 2012). In addition, the rapid increase in the global population (the population increased from 5.281 billion to 7.585 billion from 1990 to 2020, i.e., by about 1.45%) (UN, 2019) and rapid economic development (the per capita gross domestic product (GDP) increased from $4279 to $11033 during 1990-2020 and the urban sprawl increased by 1.8 times (UN, 2015)) have led to an increase in not only basic needs (e.g., food, oxygen, and water) but also a higher level of demand (e.g., dietary structure, standard of living, and spiritual fulfillment). The decline in the ES supply and the increase in the ES demand pose a significant direct threat to sustainable development. It is crucial that the global ES be mapped and the imbalance between and spatiotemporal characteristics of the supply and demand be determined to provide effective suggestions and policies for global and regional land management, urban planning, and ecological restoration.
A large number of models and methods have been developed to quantify, assess, and map the ES supply and demand in the past two decades, which can be broadly classified into the following various categories: ecological proxy modeling methods, such as the Integrated Valuation of Ecosystem Services and Trade-offs (In-VEST) models (Arcidiacono et al 2015, Kim andArnhold, 2018), value assessment methods, direct market method, public participation methods, and questionnaire survey methods. However, most of these methods are only suitable for mapping the ES supply and demand at the landscape, watershed, and even national scales due to the unavailability of data (Perennes et al 2021). It is challenging to conduct ES mapping on a global scale. Land use/ cover (LUC) is considered to be one of the vital factors influencing ESs (He et al 2019), and the different LUC types have significantly different ES capacities (Chen et al 2022). Therefore, (Burkhard et al 2012) proposed an expert knowledge-based ES matrix model, which is widely used in ES mapping studies because of its efficiency, speed, accessibility, and practicality. Based on the ES matrix model, Burkhard et al. (2012) mapped the ES supply and demand for energy, food, and water in the cities of Leipzig and Halle in central Germany from 1990 to 2007. (Jacobs et al 2015) estimated the relationship between the ESs and LUC in rapidly urbanizing regions using an ES matrix model and projected the ES changes in 2027. By comparing the LUC changes before and after a tsunami, (Kaiser et al 2013) measured the loss of ES capacity due to the tsunami using an ES matrix model. Establishing an ES matrix is a vital to step in determining the relationship between the LUC and ES capacity since it has a strong dependence on and sensitivity to the accuracy and classification of the global LUC datasets (Bian and Lu 2013). The other crucial point is that expert models struggle to deal with regional differences in ESs caused by climate and biophysical effects. Therefore, revising the ES supply and demand separately using highly relevant data is an effective means of improving the accuracy of ES mapping.
Although ES mapping can be greatly facilitated by the use of an ES model matrix (Bian and Lu 2013), it is still challenging to achieve global ES mapping with a high accuracy and multi-temporal global LUC data, which makes it necessary to implement this method using cloud platforms. The Google Earth Engine (GEE) cloud platform includes a wide range of remote sensing ecosystem detection and geographic big data cloud computing. The rapid data processing capability and the stability of the cloud platform have led to a multitude of global-scale ecological studies. Therefore, the combination of ES mapping and the GEE will be the future trend of large-scale research on ecosystems and human well-being (Farda 2017).
Based on the open and synergistic (OS) global LUC datasets from 1990 to 2020, an ES matrix model was established to map the global ES supply, demand, and relative budgets. The global net primary productivity (NPP) and enhanced vegetation index (EVI) and the nighttime light (NTL) and world population (Pop) datasets were used to correct the ES supply and demand, respectively. Finally, the spatiotemporal characteristics of the ES supply and demand in geographical sub-regions were analyzed. The goals of this study were (1) to achieve global 1×1 km resolution ES supply, demand, and relative budget mapping for 1990-2020; (2) to analyze the global spatiotemporal characteristics of the ESs and their changes of heterogeneity; and (3) to develop targeted recommendations for the management of the land and ecology in different regions.

Data processing
The selection of the global territory land cover datasets was crucial to this study. The open synergistic (OS) global land cover datasets was employed as the basis for the assessment of the ESs and was decoded using the random forest algorithm and artificial intelligence methods based on the integration of Landsat remote sense images and using a digital elevation model (DEM) as an auxiliary tool to improve the accuracy. It has a high resolution of 30 m and classifies images as cropland, forest, grassland, shrub land, wetland, water, tundra, impervious surfaces, bare land, and snow and ice. Its classification of land types is more precise than those of other global products (Zhao et al 2021). The most crucial reason for this is that the OS datasets have a long time series of 3 decades, from 1990 to 2020, and the overall classification accuracy is greater than 75%.
The ES supply capacity of the forests near the equator is normally higher than that of forests at higher latitudes due to the differences in the biophysical conditions, so further correction of the ES can significantly improve the accuracy of the mapping. Considering the availability of data, the NPP and EVI and the NTL and Pop (Marcantonio et al 2015), which are highly correlated with ESs, were selected to correct the ES supply and demand, respectively. Since global NPP, EVI, and Pop raster data are available for 2000-2020, the raster data from 1990 to 1995 were used to create the linear regression model and ordinary least squares (OLS) model. NTL datasets are available for 1995-2010, the corresponding raster data are available for 1995, 2015, and 2020, which were also used to create the linear regression model. The product names, raster resolutions, and data sources are presented in table 1. All of the datasets were resampled to a 1×1 km resolution in the GEE to reconcile the computational differences.
The research process and methodology are illustrated in figure 1. First, the OS land cover dataset with a 30 m resolution was acquired from the GEE. Second, an ES matrix model of the relationship between the LUC types and the ES capability was established. Based on the GEE, the binarization method (e.g., the forest land was set as 1 and the other land was set as 0, and then, the proportion of forest land at a 1 km 2 resolution was calculated) was used to extract and count the proportion of each type of land cover at a 1×1 km resolution and the original ES supply and demand map was obtained through superposition of the ES capacity given by the ES matrix model. Third, because the ES capacity generated by the different combinations of ecosystem types in the region had a huge discrepancy (Blasi 2013) (e.g., the ES capacity of forest-grassland-water land combination in the region is more variations from that of shrub land-water-impervious land combination as well as the different region also have significant distinction in ES capacity), the NPP and EVI and the Pop and NTL datasets were selected to correct the ES supply and demand, respectively. Finally, the spatiotemporal and imbalance characteristics of the ES capacity in the geographic sub-regions were determined, the dynamic relationship between the ESs and human activities was analyzed, and targeted recommendations were developed to guide policy development for regional ecological protection, land management, and urban planning policy in the different sub-regions. Overall, the data processing was implemented and visualized using the GEE, and partitioned statistics were applied to export the data into the local platform for further analysis.

ES matrix model
The Millennium Ecosystem Assessment (MEA) has been widely recognized and classifies the ES functions into four major categories: provisioning services, regulating services, cultural services, and support services (Norgaard 2008). However, numerous studies have shown that many of the functions of regulating and support services conflict (Zhang et al 2020, Sasaki et al 2021). Therefore, in this paper, nine regulating service indicators (e.g., global climate regulation, flood control, and pollination), eleven provisioning service indicators (e.g., agricultural production, wood production, and water purification), and two types of cultural services (e.g., recreation) were selected to represent the ES capacity.
An ES matrix between the 22 ES indicators and 10 land cover types was established, and then, the ES capacity was divided into six degrees based on expert experience: 0 -relevant capacity; 1 -low relevant capacity; 2relevant capacity; 3 -medium relevant capacity; 4 -high relevant capacity; and 5 -very high relevant capacity. The final ES capacity was obtained by weighting the proportion of land and the capacity in the ES indicators.
In this study, the method of assigning values and then correcting them to the supply and demand in the ES matrix was adopted. Based on the the ES matrix model established by interviewing hundreds of experts and the Corine Land Cover The global ES supply and demand capacities were measured using 1×1 km grid cells, and the binarization method was used to obtain the global 1×1 km resolution proportion of the land raster. Then, based on the ES matrix model assignment, the ES capacity of the grid cell was weighted and summed as follows: where E s and E d are the supply capacity and demand capacity of the ESs, respectively; n is the number of ES indicators (21 types in this study); E k is the supply (demand) capacity degree of ES type k; M k is the area of land type k; i is the number of 1×1 km grids globally; and M i is the area of grid i.

Global ES supply and demand correction
The selection of correction indicators should consider not only the high relevance to the ESs but also the availability and reliability of the datasets. Fortunately, the results of numerous studies have shown that the EVI and NPP have strong correlations with the ES supply, Pop, while the social development (replaced by Nighttime light datasets) has a robust correlation with the ES demand (Burkhard et al 2012). Related studies have discovered that the ES demand capacity and diversity are always greater in developed countries than in developing countries (Quatrini and Crossman, 2018). In addition, several studies have also shown that middleincome people are more sensitive to the ES demand (Mascayano et al 2019), and their Pop and economic development exhibit nonlinear relationships with the ES demand. The same applies to the ES supply (Zacharias and Ressi 2018). Therefore, for the correction, the Sigmod function, which is commonly used in neural networks, was adopted to converge the correction data to the range of 0.5-1.5 (its ability to influence the results was 0.5 times higher than that of the ES matrix model) through corresponding changes. In the first step, the EVI and NPP indicators and the Pop and NTL indicators were integrated as follows: where S i and Di are the integration coefficients of EVI and NPP (supply revision) and the Pop and Light (demand revision) data for raster cell i, respectively. To facilitate comparison in the time dimension,Ē,N,P, andL are the EVI, NPP, Pop, and NTL mean raster values of 3586.61, 1626.49 kg * C m −2 , 28.8 per ha −1 , and 5.473 in 2005, respectively. w 1 and w 2 are the contribution weights of the different correction coefficients. Based on a relevant previous study (Taelman et al 2016), w 1 and w 2 were set as 0.5 and 0.5 in this study, respectively. Substituting the integration factor into the changed Sigmod function converges the correction factor to the 0.5-1.5 interval: where RS i and RD i are the correction coefficients of the ES supply and demand for raster cell i, respectively. S i and D i are the integration factors. The correction coefficient CF i is converged to between 0.5 and 1.5 by adding a constant C (0.5 in this study). When D i is 1 (i.e., the Pop and NTL are the global average) in equation (5), the ES demand correction coefficient RD i is 1. As D i gradually increases, the ES correction coefficient gradually increases, and the sensitivity gradually decreases, and the maximum correction coefficient converges to 1.5. The ES supply and demand correction coefficients are multiplied by the corresponding original ES capacity to obtain the final ES capacity. The specific equations are as follows: where ST i and DT i are the ultimate ES supply and demand values of cell i, respectively. RS i and RD i are correction coefficients. T s,i and T d,i are the ES supply and demand values of cell i, respectively, which are calculated using the ES matrix models.

Spatiotemporal analysis
Based on the results of the global ES supply, demand, and relative budget mapping, the spatiotemporal and imbalance characteristics of the ES capacity in 21 geographic sub-regions (Supplementary figure. 1(available online at stacks.iop.org/ERC/4/065002/mmedia)) based on the geographic sub-region criteria in the UN proposed geographic zoning table scheme were analyzed to provide more targeted suggestions. The spatial model used to analyze the ES capacity is as follows.
(i) Center of gravity analysis. The center of gravity model can better identify the spatiotemporal characteristics of the geographic elements and has been widely used in recent years. The center of gravity model was used to analyze the movement of the center of gravity of the ES budget in each geographic sub-region. The change in the ES budget can greatly reflect the effect of human activities on the ES capacity, so the spatial mobility of human activities can be embodied by the changes in the ES budget's center of gravity. The formula is as follows: (ii) Spatial autocorrelation analysis. Spatial statistical analysis can better analyze geographically related, spatially dependent, and spatially correlated data and can reveal the spatial correlation characteristics (Mascayano et al 2019). To analyze the spatial aggregation characteristics of the ES supply, demand, and budget, a global hexagonal fishing net with a 1°basic unit was established. Moran's I local spatial autocorrelation analysis was conducted by counting the mean ES value in the fishing net in 2020. The formulas used to explore the spatial autocorrelation and imbalance of the ES capacity are as follows: where x i and x j are the ES capacity of grid cells i and j (i ≠ j), n is the number of grids, which is the global ES capacity mean. S 2 is the data sample variance, and w ij is the spatial weight matrix (constructed using the K-nearest method). I global is within the range of (−1, 1), and a positive I global indicates the existence of a certain amount of spatial aggregation. The higher the aggregation is, the closer the value is to 1. In addition, the local Moran index was utilized to analyze the local spatial heterogeneity of the ES capacity. The local Moran index identifies the H-H and L-L aggregation patterns in the local space to obtain the heterogeneity characteristics of the ES.
(iii) Clustering analysis. Four indicators (e.g., the ES supply, demand, relative budget, and rate of change) were calculated iteratively to classify the geographic sub-regions. The two most similar data points among all of the data points were combined, the similarity between the two types of indicator points was calculated, and a clustering tree was iteratively generated. This method was used to find the similarities and differences between the ES capacities of the geographic sub-regions in order to make more effective and targeted recommendations.

Results
3.1. Spatiotemporal and heterogeneity of global ES capacity 3.1.1. Spatiotemporal characteristics of the global ES supply, demand, and budget Globally, the high ES supply capacity areas are mainly located in the equatorial and 50°N regions ( figure 3), with the equatorial ES supply capacity being as high as 63.74. The ES demand capacity is distributed between 52°S and 63°N. The ES demand in the northern hemisphere is significantly higher than that in the southern hemisphere, and it is most significant within 6-50°N. The global ES budget exhibits a double peak distribution with latitude. The equatorial and 50°N regions are the ES budget surplus regions, and they have large ecological resources and play a very vital role in global climate regulation, carbon fixation, and biodiversity maintenance. In contrast, the low ES budget regions are located in the 10-40°N region, which not only includes more than 85% of the global bare land but is also the primary area of human activities and has a high ES demand.
In terms of longitude, the ES supply is more concentrated in the Western Hemisphere, with the highest ES supply of 55.43 occurring in the longitude region corresponding to the Amazon basin. However, there are two noticeable gaps in the 73-96°E region and the 14°W-15°E. The ES supply is low, and more crucially, the human activities and high demand have significantly affected the ES budget in these regions. In the past three decades, in the 10°S-60°N region, the ES budget has decreased by 75.89%, and the growth of the ES demand has been the most significant in this region.

Spatial autocorrelation characteristics of the ES capacity
The global ES supply, demand, and budget exhibit significant (P<0.05) autocorrelation and heterogeneity, with Global Moran's I values of 0.89, 0.796, and 0.816, respectively (figure 4). Spatially, the H-H aggregation areas of the ES supply are mainly distributed in the Amazon Basin, Congo Basin, Southeast Asia, and Siberia; while the L-L areas of the ES supply are mainly distributed in the Sahara Desert, central Asia, and at high northern latitudes. In contrast, the H-H regions of the ES demand are mainly distributed in Europe-western Asia, India, eastern China, and northeastern China. The global ES budget exhibits a high degree of spatial aggregation, and there are significant difference between the ES surplus regions (ESSRs) and ES deficit regions (ESDRs), indicating a high degree of imbalance between the ES supply and demand.
The degree of spatial aggregation of the ES budgets varies significantly across the geographic sub-regions. In terms of the percentages of spatially autocorrelated factors ( figure 4), Europe, Africa, and Asia generally have higher spatial imbalance characteristics, with the most pronounced imbalance occurring in southern Europe, northern Africa, and eastern Asia, which need to use ES flows more rationally to meet human benefits. The spatial aggregation of the ES budget decreased slightly during 1990-2020 in all of the geographic sub-regions (except for northern Africa), suggesting that human activities led to the spatial heterogeneity of the ES budget having a smaller impact and the decline in the imbalance being insufficiently pronounced.

ES capacities of geographic sub-regions 3.2.1. ES capacity variation statistics
The average capacity of the ES supply was 33.31 (figure 5) and the geographic sub-regions with high ES supply capacities were Melanesia, southeast Asia, the Caribbean, South America, and central Africa, with supply capacities of 75.14, 62.58, 48.27, 47.75, and 43.31, respectively. These regions are mainly distributed around the equator. The low-value regions were mainly located in parts of Africa (northern Africa (5.43), western Africa (16.68), and southern Africa (17.18)) and parts of Asia (western Asia (11.99) and central Asia (18.85)). In these regions, rainfall is scarce, the main land cover type if bare land, and ecosystems are fragile. The global ES supply increased by 0.67% per decade; however, the ES demand increased by about 3.36% per decade (figure 5). The sub-regions with high ES demands were western Europe, southern Europe, southern Asia, and eastern Asia, with capacities of 30.08, 26.21, 23.47, and 19.08, respectively. The low ES demand regions were characterized by very small populations, primarily bare land, arid climates, and individual islands far from the mainland.
The balance line between the ES supply and demand (figure 5) shows that the ES supply capacity was only less than the ES demand in southern Asia; and the ES supply exceeded the ES demand in all of the other subregions. The global ES supply was 2.405 times higher than the demand; and a significant proportion of the ES supply was located in Melanesia, central Africa, and South America. Although the ES supply exceeded the demand in most of the geographic sub-regions (average ES budget of 19.63), there were distinct differences in the ES budgets of the sub-regions, with Melanesia (69.19), southeast Asia (42.02), South America (35.47), and Central Africa (34.25) having the highest ES budget capacities. Northern Africa (1.26), western Asia (3.77), central Asia (4.90), and western Africa (4.8) had low ES budget capacities.
Influenced by the increase in the global ES demand, the ES budget of each sub-region exhibited a decreasing trend, among which the ES budget in Europe decreased the most significantly (2.11), followed by the USA (1.63), and that of Africa decreased the least (only 0.44). The global ES budget exhibited a decreasing trend and the rate of decrease increased from −1.90% to −3.12%. Only the ES budget of Melanesia increased, with a rate of −0.69%. The ES budget decreased the most significantly in southern Asia (−8.65%), western Europe (−6.61%), and central Asia (−5.25%). In contrast, the low ES budgets in northern Asia (−021%), central Africa (−0.13%), eastern Africa (−0.97%), and A.N. Zealand (−0.965%) experienced small changes due to the relatively weak human activities and low ES supply capacities of these regions.

Shifts in the center of gravity of the ES budget in the geographic sub-regions
The center of gravity coordinates, migration distance, and migration direction of the ES budget from 1990 to 2020 were calculated for each sub-region (table 2). The average migration distance of the center of gravity of the ES budget in the sub-regions was 34.17 km. It migrated 44.23 km and 36.74 km in Europe and Asia, respectively, and by only 27.518 km in Africa. The migration distances were longest in western Asia (73.5 km), central Asia (60.9 km), southern Europe (55.34 km), and western Europe (954.39 km). It was found that the sub-regions with higher rates of change in the ES demand and higher intensities of human activity perturbations generally had longer migration distances. however, the spatial distribution of the ES budget in Melanesia (4.96 km), central Africa (2.59 km), and western Africa (7.73 km) did not change significantly, indicating that the spatial perturbations of the ES capacity by human activities was not significant.
From the perspective of the migration direction, the center of gravity of the ES budget in Europe migrated to the east, while in northern Africa and central Africa, it migrated to the south and southwest, respectively. The center of gravity of the ES budget exhibited a certain pattern of migration, it tended to move away from the ESDR in the Sahara Desert and close to the ESSRs in Siberia in the European region. in northern and western Africa, it tended to migrate toward the Congo Basin ESSR, while in southern and eastern Africa, it tended to migrate toward the Congo Basin. This phenomenon also occurred in Asia. In northern and western Asia, the center of gravity of the ES budget shifted to the west and northwest, respectively, while in the other Asian sub-regions, it tended to shift to the east and southeast. It should be noted that southeast Asia is an ESSR, with abundant resources and a superior environment. The ES budgets of Central America and North America tended to move closer to the Amazon watersheds. The imbalance between the ES supply and demand produced a push-pull effect, which forced humans to move closer to the ES surplus regions and away from the ES deficit regions, and the destruction of the ecological environment promoted this phenomenon.

Human-ES dynamic and policy implications 4.1.1. Push-pull effect and human-ES dynamics
The imbalance between ES supply and demand generated a push-pull effect on humans, forcing people to move toward the ES surplus regions. To further analyze the push-pull mechanism of the ES capacity and human activities, a diagram was established to illustrate their relationship (figure 6). From the supply perspective, natural factors are the main cause of the ES supply capacity, while LULCs are a vital link between humans and ESs (Qin and Fu 2019). The push-pull effect is mainly influenced by two aspects: ecological pressure and environmental vulnerability (Dai et al 2021). The urban sprawl, forest shrinkage, and decrease in habitat quality caused by human activities produce environmental pressure. However, the ecological pressure needs to be extremely severe to force humans to move toward ES surplus regions (ESSRs) and to keep them away from ES deficit regions (ESDRs). Humans also consciously move away from ESDRs and closer to ESSRs when the ES supply (e.g., the food supply and water supply) is insufficient, (e.g., fleeing, searching for water sources, migrating to cities with better climate and environmental conditions). It should be noted that environmental vulnerability results in a region being subjected to frequent disasters such as soil erosion, landslides, floods, and others natural disasters. Environmentally vulnerable regions generally produce a push force that drives humans away from the region to obtain more convenient and high-quality natural resources and service products.
The increase in the human ES demand is mainly reflected in two crucial aspects: population increase and socio-economic development. Everyone has basic ES demands such as the need for food, water, and oxygen, which means that an increase in population will increase the basic ES demand (Daams and Sijtsma, 2013). This results in human activities being driven toward ESSRs where the ES supply is abundant and abundant resources and products can be obtained. Increasing socio-economic development leads to a diversification of the human ES demand beyond the basic ES demands, such as eating more nutritious food and drinking purer water, as well as intangible but healthier spiritual demands (e.g., forest tourism, hiking for exercise, landscape aesthetics, and history and culture). The diversity of the demand characteristics and the higher levels of demand result in the ES supply being unable to satisfy everyone's demands, forcing human activities to consciously and subjectively search for a more diverse ES supply. It is for these reasons that human activities move closer to ESSRs and farther away from ESDRs. Geographic sub-regional categorization and differential response policies Classification of the ESs in 21 geographic sub-regions was conducted through systematic clustering with the ES supply, demand, relative budget, and rate of change as the factors (Supplementary Fig. 2). It can be seen that the sub-regions can be divided into two ES macro regions, and eight ES sub-regions were classified (figure 7). Subregions I, II, III, and IV comprise one of the ES macro regions, with a better ES supply and habitat quality, while sub-regions V, VI, VII, and VIII comprise the second ES macro region, with a better ES supply and ecologically fragile and sensitive regions. Based on the specific characteristics of the eight ES sub-regions, we propose targeted recommendations to guide local governments in the development of policies for environmental management, urban planning, and ecological conservation in these regions. ES sub-region I only includes Melanesia, accounting for 0.29% of the global territory. This region is mostly composed of islands and is characterized by few human activities. Since this region is located in the equatorial region, it has the highest ES supply and budget. It can provide more ES products and plays a crucial part in atmospheric regulation, carbon sequestration, and global warming, so establishing ecological reserve areas is a vital way to prevent a decrease in the forest area.  ES sub-region II includes western and southern Europe, accounting for 2.01% of the global territory. Although these sub-regions have a high ES supply, this region also has a high ES demand and its ES budget significantly decreased during the study period. Therefore, establishing ecological reserve areas and stopping the rapid urban sprawl will ensure that the regional ES budget is in a state of surplus. The ES flow should receive more attention in these sub-regions to resolve the imbalance between the ES supply and demand.
ES sub-region III includes northern Asia, northern Europe, and Central America, accounting for 19.57% of the global territory. It has a high ES supply and a weaker ES demand, and its ES budget decreased slightly during the study period. Ecological reserve areas should be established to provide direct and indirect ES products regionally as well as globally.
ES sub-region IV includes eastern Europe, the Caribbean, South America, central Africa, and southeast Asia. These sub-regions have a high ES supply and surplus budgets, and their ES demands moderately and slightly decreased during the study period. It is also suitable to establish ecological reserve areas in these sub-regions and design the orientation of the ES product flow to cope with the imbalance between the ES supply and demand.
ES sub-region V only includes southern Asia, accounting for 3.93% of the global territory. It is an exceptional sub-region in that its ES demand exceeds its ES supply; therefore, the environment faces serious problems in this region. The quality of the environment in region V needs to improved, ecological restoration and compensation policies need to be implemented, and ecological reserve areas should be established to avoid further degradation of the environmental quality. In addition, the rapid population growth and urban sprawl should be controlled to avoid decreasing the area of cropland and forests, as well as the biodiversity.
ES sub-region VI includes western Asia, Australia and New Zealand, central Asia, western Africa, and southern Africa, accounting for 15.56% of the global territory. These sub-regions have a low ES supply and high demand. It is crucial that the ES supply is only slightly greater than the demand and that the ES budget decreased significantly during the study period. Controlling the rapid population growth and urban sprawl is an effective measure to stop the increase in the ES demand (Mendiola et al 2015). In addition, enhancing the quality of the ecological environment, implementing ecological restoration and compensation policies, and sustainably managing ecological resources are effective measures of attaining human benefits.
ES sub-region VI includes eastern Asia, North America, and eastern Africa, accounting for 33.98% of the global territory. It has a moderate ES supply and demand and the rate of change of its budget is low. Although the ES supply meets the ES demand in these sub-regions, theoretically, the spatial heterogeneity of the ESs forces us to consider how to control and manage the flow of the ESs. Establishing ecological reserve areas will help achieve sustainable development, and ecological restoration and compensation policies should be adopted.
ES sub-region VIII only includes northern Africa, accounting for 3.42% of the global territory, and the LULC types in this sub-region are mainly bare land and the Sahara Desert. Therefore, this region is characterized by drought, ecological vulnerability, and a low ES supply capacity. Considering the small amount of human activity, afforestation, and soil erosion prevention policies should be adopted in these regions, and a core ecological protection zone should be established to prevent the expansion of desertification.

Ecosystem service budget mapping
We discovered significant spatial heterogeneity in global ecosystem service (ES) budgets (Global Moran's I value is 0.816), which greatly intensifies the inequity of development and human well-being. Although the conclusion that ES supply is 2.405 times greater than demand, which is an optimistic figure, concentrated ecosystem service deficit regions such as the Sahara desert to Mongolia zones impede ecosystem flows.
Land cover data has been widely used to map ES at a local scale, which supports decision-making for ecological conservation and management. As for global ES mapping, cloud computing support, such as the Google Earth Engine cloud platform, is needed to improve computational efficiency. Therefore, this work can be a very meaningful and successful attempt. In addition, the obtained global ES budget maps with high spatial resolution in our study provided an opportunity to promote the strengthening of ecological conservation cooperation among various countries, enhance the fairness of stakeholders, optimize government ecological management decisions, and ultimately achieve the sustainability of the environment.

Limitations and uncertainties
An ES matrix model based on expert empirical knowledge was adopted to map the global relative ES budget. Therefore, the ES model has a strong sensitivity to the classification accuracy and classification type the LUC datasets, and the LUC classification has the greatest impact on the accuracy of the ES budget. (Burkhard et al 2012) scored the ES capacity of 44 different types of the coordination of information on the environment (CORINE) land use/cover data. It is not surprising that the use of coarse global LUC data decreases the accuracy of the relative ES budget and causes limitations since more land types allow consideration of the subtle changes in carbon sequestration and water purification and also take the combination characteristics among the different land types into account.
The expert-based matrix model is used to reflect the relative ES budget capacities through assignment of the ES capacity of the land, so it is only a semi-quantitative model (Burkhard et al 2012). This method is highly subjective, operational, and practical; however, it needs to be improved in terms of credibility, consistency, and validity. All of the ES budget results are relative values, but it can still be concluded that the vast majority of the geographic sub-regions maintained a state in which the ES supply is greater than the ES demand during the study period (except for southern Asia, which exhibited the opposite phenomenon). Because the ES matrix model only defines cropland land and impervious surfaces as ES deficit land and on the sub-region scale the ES deficit areas (impervious surfaces and cropland) accounted for a small percentage, the results need to be verified further in future studies. The use of the EVI, NPP, and Pop datasets and the NTL dataset to correct the ES supply and ES demand has limitations and is uncertain because the ES capacity is influenced by a variety of factors, such as socio-cultural distinctions, the age and gender structure of the population, and human perceptions of the environment. Regarding the push-pull effect, only the effect of the environment on human migration was considered. Although socio-economics are also a vital factor, they were not taken into consideration in this study. Furthermore, due to the unavailability of datasets, it was difficult to consider all of the relevant factors when correcting the ES capacity. Therefore, representative factors were employed to improve the accuracy of the relative ES budgets. it should be noted that the global ES supply was calculated to be 2.405 times the ES demand, which differs from Millennium Ecosystem Assessment (MEA) for multiple reasons; however, the imbalance between the ES supply and demand and the relative budgets also provide a scientific basis for developing suggestions globally.

Conclusions
Based on open and synergistic (OS) land use/cover datasets and the GEE, an expert-based ES matrix model was employed in this study to map the global ES supply, demand, and relative budget from 1990 to 2020. The results of this study revealed that the global ESs exhibited strong heterogeneity and the ES supply was 2.405 times higher than the demand. We also discovered that the growth rate of the ES demand was 3.36% per decade. Western and southern Europe and southern and eastern Asia were the high demand capacity areas, among which southern Asia's ES demand exceeded its ES supply. The imbalance between the ES surplus regions and ES deficit regions produced a push-pull effect, which forced humans to move toward ESSRs and away from the ESDRs. Moreover, the destruction of the ecological environment will accelerate this phenomenon. Eight ES sub-regions were defined, and targeted land management, urban planning, and environmental remediation policies should be established for each of these regions.