*2.1. Study Area*

The study area was in the middle and lower reaches of the Yellow River basin (103◦36 – 119◦55 E, 41◦03 –32◦46 N), at an altitude of about 0–4082 m. The area was located in the central-eastern part of China (Figure 1). The main provinces involved included Inner Mongolia Autonomous Region, Henan Province, Shaanxi Province, Shanxi Province, Qinghai Province, and Gansu Province. The middle and lower reaches of the Yellow River Basin are dominated by plains and hills. The region has a temperate continental climate and a temperate monsoon climate. The middle and lower reaches of the Yellow River Basin have abundant light, high temperature, abundant precipitation, are suitable for crop growth, and are the main production areas for agricultural products. The Yellow River Basin is an important economic zone and an important base for energy, chemicals, raw materials, and basic industries in China.

#### *2.2. Data Sources*

In this study, we used data from five different sources. (1) Meteorological elements and daily precipitation for 2000, 2010, and 2018, supplied by the China Meteorological Data Sharing Network (http://data.cma.cn/ (accessed on 1 July 2021)), were batch interpolated using the professional meteorological interpolation software ANUSPLIN. (2) Monthly NDVI data for 2000, 2010, and 2018 at a spatial resolution of 1 km × 1 km, supplied by the Geospatial Data Cloud (http://www.gscloud.cn/ (accessed on 1 July 2021)), and annual NDVI data were obtained using the maximum synthesis method. (3) Population data for 2000, 2010, and 2018 at a spatial resolution of 1 km × 1 km was supplied by WorldPop (https://www.worldpop.org/ (accessed on 1 July 2021)). (4) Land use data for 2000, 2010, and 2018 at a spatial resolution of 1 km × 1 km was supplied by the Resource and Environmental Science and Data Centre (http://www.resdc.cn/ (accessed on 1 July 2021)). (5) Grain production, energy consumption, and water consumption by the municipality for 2000, 2010, and 2018 was obtained from the statistical yearbooks and water resources bulletins of each province.

**Figure 1.** (**a**–**c**) Location, elevation, and 2018 land use of the study area.

#### *2.3. Framework for Optimizing Production-Living-Ecological Space*

Based on previous approaches and frameworks for land use optimization [23,44], this study identified ideal land-use patterns and optimized PLES at different scales by quantifying the mismatch between the supply and demand of ES associated with PLES. This was achieved through the following four core steps (Figure 2): In the first step, the composition, configuration, and spatial transition of PLES were analyzed based on land use data in the Yellow River Basin during 2000, 2010, and 2018. This step aims to examine the spatial changes in PLES in the Yellow River Basin during 2000, 2010, and 2018, and establish a basis for subsequent research. In the second step, those ES suitable for the Yellow River Basin were selected based on the basic principles for the selection of ES proposed above to assess the ES supply and demand. Mismatches and shortages between ES supply and demand were also identified. In the third step, based on the correlation between the ratio of production/living/ecological space and the supply and demand of ES, the thresholds were identified when ES supply and demand were imbalanced. This step aims to analyze the links that exist between the two, the most central part of which is the identification of thresholds. In the fourth step, the thresholds identified in the previous step were used to identify optimization areas using ternary phase diagrams, which were then optimized for PLES at different scales.

Step 1: The classification of production-living-ecological space based on land use types.

Production space is mainly the area that provides various products or services for people. Living space refers to the area that provides the function of carrying and guaranteeing human habitation and provides the function of residence, consumption, leisure, and recreation in the country. Ecological space refers to the area that can provide an ecological barrier and has the function of regulating the atmosphere, concealing water, and maintaining soil and water [11,12]. In this study, a classification system for PLES in the Yellow River

Basin was constructed based on geographical features and previous research results [45] (Table 1).

**Table 1.** Production-living-ecological space classification in the Yellow River Basin.


**Figure 2.** Optimization framework for production-living-ecological space.

Step 2: Quantify ES supply and demand and identify spatial mismatches.

(1) Quantifying ES supply and demand

Water yield service are the ability of an ecosystem to intercept or store water resources from rainfall while mitigating ground runoff [46]. The Yellow River Basin is an important water source in northwest China and assessing the water yield service is of great practical importance for the rational use and conservation of water resources. The water balance model was used to calculate the supply of water yield service [47]. The amount of water consumed per capita in each city in the Yellow River Basin was obtained from data on industrial, agricultural, and domestic water consumption and the resident population of each city and combined with data on population density to obtain the demand for water yield services. The formulas are as follows.

$$WY\_{(x)} = P\_{(x)} - ET\_{(x)} \tag{1}$$

$$ET\_{\langle x\rangle} = \frac{P\_{\langle x\rangle} \left(1 + PET\_{\langle x\rangle} / P\_{\langle x\rangle}\right)}{1 + \omega\_{\langle x\rangle} \left(PET\_{\langle x\rangle} / P\_{\langle x\rangle}\right) + \left(PET\_{\langle x\rangle} / P\_{\langle x\rangle}\right)^{-1}} \tag{2}$$

$$D\_{\rm w} = D\_{\rm pw} \times \rho\_{\rm pop} \tag{3}$$
 
$$\text{quad } \text{uvorravo vold } \text{at } \text{nivol } \times (\text{m}^3); D\_{\rm v} \text{ : ir } \text{the } \text{nivorravo vold}$$

where *WY*(*x*) represents the annual average yield at pixel x (m3); *P*(*x*) is the average annual precipitation on pixel x (mm); *ET*(*x*) is the actual annual evapotranspiration at pixel x (mm); *PET*(*x*) is based on the Penman-Monteith formula [48]. *Dwy* is water demand, which in this case equates to water consumption (m3); *Dpw* is water consumption per capita (m3/ person), which includes water consumption for agricultural, industrial, domestic, and ecological purposes; *ρpop* is the local resident population density (person/km2).

Grain production service, as one of the basic ecological services, plays a vital role in human survival and development [49]. There is a significant linear relationship between crop and livestock production based on the NDVI. The total production of grain was allocated according to the ratio of raster NDVI values to total arable land NDVI values, which in turn characterized the grain production capacity of each raster. Grain demand was estimated by multiplying the per capita grain demand by the population density [17]. The formulas are as follows:

$$GP\_{(x)} = GP\_{sum} \times \frac{NDVI\_{x}}{NDVI\_{sum}} \tag{4}$$

$$D\_{\mathcal{S}} = D\_{\mathcal{PS}} \times \rho\_{\text{pop}} \tag{5}$$

where *GP*(*x*) is the total production of grain for grid x (t/km2); *GPsum* is the production of grain products for each province (t); *NDV Ix* is the normalized difference vegetation index for grid x; *NDV Isum* is the sum of NDVI of cropland for each province; *Dg* is the grain demand (t/km2); *Dpg* is the annual per capita grain consumption (t/person); *ρpop* is the resident population density (person/km2).

Carbon sequestration services are important regulatory services in ecosystems. The CASA model is a common model for calculating NPP due to its high calculation accuracy and easy-to-access data and parameters [50]. The carbon sequestration demand was calculated from the product of population density, per capita energy consumption, and energy carbon conversion rate, where energy consumption was obtained from the statistical yearbooks of the Yellow River Basin provinces. The formulas are as follows:

$$NPP(\mathbf{x},t) = APAR(\mathbf{x},t) \times \varepsilon(\mathbf{x},t) \tag{6}$$

$$D\_{\mathfrak{c}} = D\_{\mathfrak{pc}} \times \rho\_{\text{pop}} \tag{7}$$

where *NPP* is the net primary productivity of the pixel x at time t (gC/m2·a); *APAR* is the Absorbed Photosynthetic Active Radiation (MJ/m2·a), which is estimated from the ratio of total solar radiation (SOL) to absorbed photosynthetically active radiation (FPAR); *ε* is the efficiency of conversion of photosynthetically active radiation to organic carbon (gC/MJ2), which is estimated by maximum light energy utilization (0.389 gC/MJ2), temperature stress (Tε), and water stress (Wε). *Dc* represents carbon sequestration demand (t/km2); *Dpc* is the annual per capita carbon consumption (t/person); *ρpop* is the resident population density (person/km2).

Soil conservation service reduces soil erosion and restores soil fertility, which is critical to agricultural production [51]. This study quantified the soil conservation service supply based on the classical revised universal soil loss equation (RUSLE). The ecosystem service demand is the number of ecological goods that humans expect to be able to obtain from an ecosystem. Since actual soil erosion causes unwanted human losses and humans expect to manage these actual amounts of soil erosion, the actual amount of soil erosion is defined as the soil conservation service demand. The formulas are as follows:

$$SC = R \times K \times LS \times (1 - C \times P) \tag{8}$$

$$D\_s = R \times K \times LS \times \mathbb{C} \times P \tag{9}$$

where *SC* is soil conservation; *Ds* is actual soil erosion; *R* is the precipitation erosion factor; *K* is the soil erosion factor; *LS* is slope length factor; *P* is soil conservation factor; *C* is vegetation cover factor [52].

(2) ES supply and demand mismatches and shortfalls

The supply and demand of ES are significantly spatially heterogeneous and are reflected in spatial mismatches. The state of ES supply and demand can be characterized by the ecological supply-demand ratio (ESDR), which can be used to reveal the nature of surpluses or deficits [35,53].

$$ESR = \frac{ESS - ESD}{(ESS\_{max} + ESD\_{max})/2} \tag{10}$$

where *ESS* and *ESD* refer to the ES supply and demand, respectively; *ESSmax* is the maximum value of the ES supply; *ESDmax* is the maximum value of the ES demand. *ESDR* > 0 indicates a surplus, *ESDR* = 0 indicates a balanced ES supply and demand, and *ESDR* < 0 indicates a deficit.

Step 3: The impact of production-life-ecological space on the ES supply and demand imbalance.

The ESDRs for the four major ES in 2000, 2008, and 2018 were calculated by the above method, while the production space ratio/living space ratio/ecological space ratio at the 1 km grid scale was calculated based on 30 m land use data. The data were statistically graded for the years 2000, 2010, and 2018, and then least squares regression analysis was conducted via SPSS to plot the trend line between ESDR and production space ratio/living space ratio/ecological space ratio to indicate negative or positive effects and significance levels. Spatial land management thresholds (i.e., the ratio of production-living-ecological land when there is a deficit in the ES) were then calculated based on the results of the regression analysis.

Step 4: Identification of the direction of optimization and policy recommendations.

A ternary diagram is a type of center of gravity diagram that has three variables but requires the sum of the three to be constant. In an equilateral triangular coordinate system, the position of a point in the diagram represents the proportional relationship between the three variables. In this study, the same ternary was used to visually express the ratio of production-living-ecological space, which was used to identify the optimization area with the main optimization direction, where the ratio occupied by this type of land at the endpoint is 100%. The regions are divided according to the thresholds determined above. The projection of units of different scales is performed, and when the projection falls in the ideal region, it means that the unit does not need to be optimized, and when it falls in other regions, the direction and quantity relationship of optimization can be determined based on the direction and distance from the ideal region.

#### **3. Results**
