Investigating the Coupling of Supply and Demand for Urban Blue and Green Spaces’ Cooling Effects in Shandong, China

: It is of great signiﬁcance to determine the level of demand for thermal environment regulation and the availability of blue–green spaces for thermal environment regulation to alleviate the effects of urban heat islands. Taking Shandong Province, China, as the study area, combined multi– source remote sensing data are used in this study to construct the index system of cold island supply capacity (CIS) and the cold island demand level (CID). We use the methods of spatial regression, quadrant division, and coupling coordination degree to analyze the correlation, matching status, and the level of coordinated development between the supply capacity and demand for the cooling effect. We also explore the change law and spatial characteristics of the blue–green spaces’ cooling effects supply and demand matching. Results show that: (1) The CIS and the CID are signiﬁcantly negatively correlated and spatially heterogeneous in distribution, with a signiﬁcant spatial spillover effect. (2) The dominant type of supply and demand is one of low supply and high demand, which means that the supply and demand for cool islands’ cooling effect are unbalanced, with signiﬁcant problems of spatial mismatch and quantitative imbalance. (3) The coupling between supply capacity and demand level is on the verge of becoming dysfunctional because the uneven distribution of urban buildings, population, and blue–green spaces reduce the coupling between supply and demand levels. This research can provide a new perspective and scientiﬁc basis for the study of the cooling effects of blue and green spaces and the mitigation of the heat island effect in densely populated urban centers.


Introduction
China's urbanization process as a whole is on a rapid upward trend.According to the latest forecast, China's urban population will reach 1.12 billion by 2050 [1].Continued urban population growth drives rapid urban expansion, resulting in the replacement of natural surfaces such as vegetation and water bodies with impervious water and the continuous reduction in blue and green spaces with natural cooling effects [2].Rapid changes in the urban landscape and increased emissions from human-made heat sources are leading to increasingly pronounced urban heat island effects [3,4], which have not only destroyed ecological balance, increased energy consumption, and environmental pollution but also pose a significant threat to human health [5,6].Studies have proven that urban heat can trigger a range of physical and psychological disorders [7].During heat waves, the incidence of respiratory and cardiovascular diseases increases significantly, and the health risks for urban dwellers continue to rise [8].Therefore, how to mitigate the adverse effects of heat islands and improve the thermal comfort of residents has become a major research focus.
Researchers have developed a range of cooling strategies for managing urban heat islands [9].Part of the research focuses on mitigating solar radiation using innovative reflective materials for pavements, roofs, and facades.In addition, shading devices and water-based artificial facilities are built into the urban design to address the effects of heat islands.There is also research into the construction of urban ventilation corridors to reduce urban heat.However, these solutions inevitably face limitations [10] related to residents' unfamiliarity with innovative building materials and techniques, the effectiveness of reflective materials that may be influenced by local climate, and the similarity of the cooling effect of urban ventilation corridors to the wind environment [11].
Impervious surface areas and blue-green spaces are important urban land use/cover patterns and structure characteristics.They can have different, even opposite, effects on how the surface radiation and energy are distributed, thus having a significant impact on the development and effects of urban heat islands [12,13].The high-temperature heat island areas and the areas with high impervious surface coverage have significant consistency in spatial distribution [14], and the more complex the patch boundary shape is, the higher the heat island intensity is [15].Impervious surface changes the reflectivity and roughness of the ground surface, affecting the heat exchange and turbulence transport between the ground and the atmosphere [16,17], altering the local microclimate and, thus, exacerbating the heat island effect.The difference in thermal properties of blue and green spaces caused by the difference in their underlying surface can effectively mitigate the thermal environment through their cooling effect [18].Green spaces absorb heat through shade and transpiration [19,20], and the water bodies have a significant cooling capacity due to their large specific heat capacity [21,22].Scholars have extensively studied the relationship between blue and green spaces and surface temperature.They found that the cooling effect closely relates to its area, shape, landscape pattern, and vegetation type [23][24][25][26].For example, the average patch area and maximum patch index of green spaces are significantly negatively correlated with its ambient temperature [27], but the patch size is not linearly related to the cooling effect, suggesting there may be an optimal patch area [28].The green space coverage is negatively correlated with temperature.When its coverage exceeds 30%, it has a noticeable mitigation effect on the heat island.When the coverage exceeds 50%, the cooling effect is considerable [29,30].The shape of green space is negatively correlated with the cooling effect [31], and the cooling effects of different green space types are also somewhat different (arbors > shrubs > grasslands) [32].The cooling effects of the water bodies have similar characteristics to green spaces.The larger the proportion of the water surface is, the higher the cooling intensity is; the larger the shape index is, the wider the cooling range is.The cooling effect of the centralized layout of a continuous body of water is stronger than that of the decentralized layout [33].Meanwhile, the environment around the water body also influences its cooling effect, which gradually weakens with the increase in density of the surrounding buildings [34].In addition, green spaces and water bodies have a synergistic effect in reducing surface temperature.This is especially manifested when the coordinated organization of blue and green spaces exceeds the superimposed cooling effect of a single green area or water body [35].Therefore, the rational organization of the morphological structure of green spaces and water bodies is a crucial factor influencing heat island patterns [36].
The current research mainly focuses on the mechanisms of land cover change, bluegreen landscape patterns, and urban spatial structure influencing the urban thermal environment [37][38][39][40].Although much research progress has been made, a number of issues still need improvement.First, research exploring urban thermal regulation services from the perspective of supply and demand is not rich enough [41].Second, the supply and demand relationship between the blue-green spaces' cooling effect and thermal regulation services is unclear [42].Third, the selection of indicators focuses on biophysical factors, while socioeconomic and other human factors closely related to heat islands are less considered.Therefore, taking China's Shandong province as the study area, we constructed the evaluation index system to quantify the cold island supply capacity (CIS) and cold island demand level (CID).The specific purposes of this study are (1) to explore the spatial relationship and distribution characteristics of cooling effect supply capacity and demand level in thermal environment regulation service, (2) to clarify the supply and demand relationship and the level of coupling in the coordinated development of the cooling effect of blue-green space, and (3) identify the key areas of imbalance between supply and demand to arrive at the accurate configuration of the artificial, built environment, and blue-green spaces.This research provides a new perspective and direction for urban development, construction, and mitigation of heat island effects.

Study Area
Shandong Province is located in the lower reaches of the Yellow River, bordering the Bohai Sea and the Yellow Sea in the east, with a land area of 155,800 km 2 and a population of nearly 102 million in 2020 (Figure 1).The province has a monsoon climate of medium latitudes with an average summer temperature of 24 • C-28 • C. The terrain is mainly mountainous and hilly, with the Shandong Peninsula in the east, the North China Plain in the west and north, and mountains and hills in the south and center, forming a landscape with the mountains and hills as the backbone and the plains and basins interspersed among them.The unique spatial patterns of the ground influence the characteristics of the urban thermal environments.According to the statistical yearbook of Shandong Province, in 2020, its population was the second largest in China, the area of urban construction land was the second largest in China, the green space coverage was the 15th largest in China, per capita green space was only 17.7 m 2 , and the total energy consumption was the first in China.Because of the rapid urban expansion, increased population concentration, and limited space for construction planning and development, the heat island effect has been intensifying, becoming a significant constraint in achieving green and high-quality development in Shandong province.
velopment, construction, and mitigation of heat island effects

Study Area
Shandong Province is located in the lower reaches of the Bohai Sea and the Yellow Sea in the east, with a land area of 15 of nearly 102 million in 2020(Figure 1).The province has a m latitudes with an average summer temperature of 24 °C-2 mountainous and hilly, with the Shandong Peninsula in the in the west and north, and mountains and hills in the south scape with the mountains and hills as the backbone and the pl among them.The unique spatial patterns of the ground influe urban thermal environments.According to the statistical year in 2020, its population was the second largest in China, the are was the second largest in China, the green space coverage w per capita green space was only 17.7 m 2 , and the total energy in China.Because of the rapid urban expansion, increased po limited space for construction planning and development, th intensifying, becoming a significant constraint in achieving gr opment in Shandong province.

Datasets
Landsat8 OLI/TIRS remote sensing image data from (http://www.gscloud.cn/(accessed on 12 October 2022)), ima with a spatial resolution of 30 m, a revisit period of 16 days, 10% were used in the study.See Appendix A (Table A1) for mospheric correction method based on a radiative transfer mo

Datasets
Landsat8 OLI/TIRS remote sensing image data from the Geospatial Data Cloud (http://www.gscloud.cn/(accessed on 12 October 2022)), imaged in the summer of 2021, with a spatial resolution of 30 m, a revisit period of 16 days, and cloud cover of less than 10% were used in the study.See Appendix A (Table A1) for image information.The atmospheric correction method based on a radiative transfer model was used, and Landsat8 TIRS 10band was selected to retrieve the land surface temperature (LST) [43].Population density and GDP data were downloaded from the Resource and Environment Science and Data Center (https://www.resdc.cn/(accessed on 12 October 2022)), 2019, with a resolution of 1 km.ZY-3 satellite imagery was taken from the general monitoring station with a resolution of 2 m for panchromatic images and 6 m for multispectral images.Night light data were from the NPP satellite VIIRS sensor with a resolution of 500 m, and POI data were from Baidu Map API interface, crawled through Python.We used night lighting and POI data to extract the boundaries of urban areas.At the same time, ArcGIS sky maps and urban and rural construction statistical yearbooks were used to optimize built-up area boundaries.ZY-3 satellite imagery combined with ArcGIS maps was used to identify blue-green spaces within urban built-up areas through automated computer extraction and manual visual interpretation.We used ENVI5.3software and ArcGIS10.2software to carry out preprocessing work, such as geometric correction, radiometric calibration, atmospheric correction, and image cropping on multisource remote sensing data.

Ecological Livability Assessment
Ecological livability is an index that characterizes an ecological function [44], represented by two characteristic indicators of green space area and distribution.The green space rate index mainly reflects the ratio of the total area of woodland, grassland, and other types of green space in the built-up areas of the city to the entire built-up area, which can indirectly point to the potential resource capacity of public green space that residents can enjoy.Green park space accessibility index refers to the ratio of the area covered by the ten-minute walkable range around the park square green space in the built-up areas of the city to the entire built-up areas, which can indicate the convenience and universality of public access to public green space, with the following formulas: where UGR is the green space ratio index of the built-up areas, A UGR is the normalization coefficient of the green area ratio index of the built-up areas with the reference value of 182.4125,UGRA is the total area of all types of green spaces in built-up areas, UA is the total of the built-up areas, UPR is the green space accessibility index of the built-up areas, A UPR is the normalization coefficient of the green space accessibility index of the built-up areas with the reference value of 111.1111, and UGA is the 800 m range around the park square green space.

LST Retrieval
This research uses the atmospheric correction method based on the radiative transfer model to retrieve the surface temperature.The radiative transfer equation is first established to obtain the radiation intensity obtained by the sensor, from which the brightness temperature corresponding to the thermal radiation intensity is calculated and finally converted to the actual temperature of the ground surface.The expression for the brightness value L λ of the thermal infrared radiation received by the satellite sensor is where ε is the surface-specific emissivity, T S is the real temperature of the ground surface K, B(T S ) is the blackbody radiation brightness W/(m 2 ×µm× sr), τ is the transmission rate of the atmosphere in the thermal infrared band, L ↑ is the upward radiation brightness of the atmosphere, and L ↓ is the downward radiative brightness of the atmosphere.According to Equation (4), the blackbody radiation brightness (Equation ( 5)) can be obtained: T S can be obtained using Planck's formula, expressed by Equation (6) as For TIRS Band10, K1 = 774.89W/(m 2 ×um×sr), K2 = 1321.08K.The above algorithm requires two parameters: surface-specific emissivity and atmospheric profile parameter, which are calculated in ENVI using the band math tool.

Cooling Effect Supply Capacity and Demand Level Index System
The CIS reflects the cold island supply capacity to mitigate the urban thermal environment and is mainly related to the cold island intensity, blue-green infrastructure, and economic development.The research shows that the cold island patch aggregation can provide a good cooling effect.Specifically, with the decrease in LST of the cold island patches, the supply capacity of the CIS is greater [45].Meanwhile, the cooling effect is related to the area and shape of green spaces and water bodies, so the larger the area is, the higher the cooling benefit is, and the more complex the shape is, the lower the heat island intensity is [46].The blue-green space is a critical vehicle for urban heat island mitigation through shading, transpiration, and heat absorption, which create the "cold island effect" [47].In addition, GDP is also a key indicator, with higher GDP providing more financial resources to mitigate the heat island effect [48].Therefore, the CIS comprehensive evaluation index system consists of the four main indicators shown in Table 1.The CID reflects the degree of cold island demand to mitigate the urban thermal environment, mainly related to the heat island intensity, population density, and the proportion and shape of impervious surfaces.Heat island patch aggregation is the primary source of urban thermal environmental risk, characterizing physical features such as risk presence and distribution [49].Population density reflects the key actors of the urban thermal environment-social attributes such as the number and aggregation degree of people exposed to thermal risks [50].Replacing a vegetated water body with an impervious surface changes the surface reflectivity and roughness, affecting heat turbulence transport and exchange, thus leading to increased demand for cooling effects.Therefore, the comprehensive CID evaluation index comprises four main indicators shown in Table 2.

Indicator Extraction and Weight Processing
This work uses ArcGIS 10.2 to create a 150 m × 150 m grid and uses grid cells to extract the evaluation values of all indicators to realize data gridding.To eliminate the dimensional and order-of-magnitude differences among various indicators, the min-max normalization method was used to normalize the positive and negative indicators in the CIS and the CID indicator system: where y + and y − are the normalization values of the positive and negative indicators in the original index value, and x max and x min are the maximum and minimum values of the original index.In this research, hierarchical analysis and entropy value methods are used to determine the weights of each indicator in the CIS and CID index systems.First, the relative importance of these measures is determined according to the decision-making experience, and a judgment matrix is formed to calculate the subjective weights.Second, the entropy method is used to calculate the objective weights.Third, the combined weights are calculated for each evaluation indicator.Finally, the CIS evaluation index system and CID evaluation index system are established: where w i and x i are the weights and normalization values of the i th indicator in the CIS index system, w j and x j are the weights and normalization values of the j th indicator in the CID index system, m is the number of indicators in the CIS, and n is the number of indicators in the CID.

Analysis of the Spatial Relationship between the CIS and the CID
The Z-score standardization method is used to match the supply and demand of the CIS and the CID and divide them into quadrants, where the x-axis represents the supply capacity of the cold island after standardization, and the y-axis represents the demand level of the cold island after standardization.The two form the four quadrants, where the first quadrant is the high supply-high demand state, the second quadrant is the low supply-high demand state, the third quadrant is the low supply-low demand state, and the fourth quadrant is the high supply-low demand state.The Z-score standardization formula is where x is the supply capacity and demand capacity of the cold island after standardization of the study unit, x i is the cold island supply capacity and demand capacity of the i th study unit, x is the mean value of the study area, and s is the standard deviation of the study area.
The spatial regression model is used to investigate the correlation between the CIS and the CID.It is expressed by Equation ( 12) as where y represents the matrix of LST dependent variables, x represents the matrix of independent variables, y represents the matrix of dependent variables, β represents the parameter matrix, λ represents the regression coefficients of the spatial residual terms, w u represents the spatial neighborhood weight matrix, and ε represents the vector of spatial error terms.
The bivariate global Moran's I is used to quantitatively characterize the spatial distribution characteristics and spatial dependence of the CIS and the CID.It is expressed by Equation ( 13) as where n represents the number of image elements in the study area, x i and y j represent the normalization values of the i th and j th CIS and CID, respectively, x and y are the mean values of the CIS and the CID in the study area, and w ij is the spatial weight matrix.
The coupled coordination model is used to analyze the level of coordinated development.The degree of coupling refers to the dynamic correlation between two or more systems that interact and influence each other to achieve coordinated development and can reflect the degree of interdependence and mutual constraints between systems.The degree of coordination refers to the degree of benign coupling in relationship interaction, which reflects the state of coordination.The coupling coordination model involves the calculation of a total of three index values, namely the coupling degree C value, the coordination index T value, and the coupling coordination degree D value.The higher C value indicates a stronger degree of interdependence and interconnection, the higher T value indicates better coordination, and the higher D value indicates better coupling and coordination between CIS and CID [51].The formula is as follows, when n = 2, assuming max U i is U 2 : when n is the number of subsystems, U i and α i are the normalization values and weights of the i th subsystem, C is the degree of coupling, and D is the degree of coordinated development.In this study, the CIS and the CID are considered equally important, and therefore, both α 1 and α 2 are taken as 0.5.According to existing research [52], the D can be divided into ten types: extreme incoordination (0.0 < D ≤ 0.1), severe incoordination (0.1 < D ≤ 0.2), moderate incoordination (0.2 < D ≤ 0.3), mild incoordination (0.3 < D ≤ 0.4), borderline incoordination (0.4 < D ≤ 0.5), barely coordination (0.5 < D ≤ 0.6), primary coordination (0.6 < D ≤ 0.7), moderate coordination (0.7 < D ≤ 0.8), good coordination (0.8 < D ≤ 0.9), and quality coordination (0.9 < D ≤ 1.0).

Ecological Livability Assessment
On the basis of the urban built-up area boundary and blue-green space data extraction in Shandong Province, the green space rate index, green park space accessibility index, and ecological livability index were calculated using Equations ( 1)-( 3).As can be seen from Figure 2, Binzhou has the highest green space index, followed by Jinan and Rizhao, while Heze has the lowest, only one-half of Binzhou's.Rizhao has the highest green park space accessibility index, followed by Dongying, while Heze has the lowest.It is worth noting that their green park space accessibility index is significantly lower than that of Dezhou, Linyi, and Liaocheng, indicating that the spatial layout of green space patches in their built-up areas is relatively fragmented and scattered, lacking a perfect blue-green space network system.The accessibility of blue-green space resource allocation and layout based on population demand is low.Binzhou has the highest ecological livability index, and Heze has the lowest.The ecological livability is divided into three grades, among which the cities with good ecological livability are Rizhao, Binzhou, Jinan, Zibo, and Dongying, the cities with medium livability are Zaozhuang, Liaocheng, Dezhou, Qingdao, and Jining, and the cities with poor livability are Weihai, Yantai, Linyi, Weifang, and Heze.On the basis of the above analysis, in order to reduce redundant work, Jinan, Qingdao, and Heze are selected as representatives.To further explore the supply-demand relationship of the cooling effect of the blue-green spaces in the thermal environment regulation service, we studied the spatial distribution patterns of the thermal environment.

Ecological Livability Assessment
On the basis of the urban built-up area boundary and blue-green space data extraction in Shandong Province, the green space rate index, green park space accessibility index, and ecological livability index were calculated using Equations ( 1)-( 3).As can be seen from Figure 2, Binzhou has the highest green space index, followed by Jinan and Rizhao, while Heze has the lowest, only one-half of Binzhou's.Rizhao has the highest green park space accessibility index, followed by Dongying, while Heze has the lowest.It is worth noting that their green park space accessibility index is significantly lower than that of Dezhou, Linyi, and Liaocheng, indicating that the spatial layout of green space patches in their built-up areas is relatively fragmented and scattered, lacking a perfect blue-green space network system.The accessibility of blue-green space resource allocation and layout based on population demand is low.Binzhou has the highest ecological livability index, and Heze has the lowest.The ecological livability is divided into three grades, among which the cities with good ecological livability are Rizhao, Binzhou, Jinan, Zibo, and Dongying, the cities with medium livability are Zaozhuang, Liaocheng, Dezhou, Qingdao, and Jining, and the cities with poor livability are Weihai, Yantai, Linyi, Weifang, and Heze.On the basis of the above analysis, in order to reduce redundant work, Jinan, Qingdao, and Heze are selected as representatives.To further explore the supply-demand relationship of the cooling effect of the blue-green spaces in the thermal environment regulation service, we studied the spatial distribution patterns of the thermal environment.

Spatial Distribution Patterns of LST
The heat island intensity is divided into five grades by mean standard deviation method.It uses the combination of average temperature and multiples of different standard deviations to define an urban heat island, which can accurately reflect the change in heat island intensity and avoid the difference in different phases [53].Figures 3 and 4 show that the portion of heat island patches in Jinan is 29.92%, mainly distributed in the central western and northeastern areas of the capital city, and the portion of cold island patches is 28.47%, distributed primarily in Daming Lake Park and Huashan Scenic Area.The hot island patches take up 28.68% of Qingdao, distributed in the south of Jimo District and the east of Chengyang District, while cold island patches take 28.46%, mainly distributed in the ecological, scenic areas such as Fufeng Mountain and Shimei Temple Park.The share of hot island patches in Heze is the highest, taking 31.19% of the city, distributed in the central urban area and the west of Mudan District, while the proportion of cold island patches is 29.61%, distributed in the northeast and southeast regions of the urban area.In terms of spatial distribution, the heat island areas are concentrated in the central urban areas where the proportion of impervious surface is high and the vegetation coverage is low, especially in the industrial and commercially packed areas with high population density, complex built environments, and scattered the green space patches.However, cold islands are mainly located in the suburbs and outer suburbs.The suburban area is adjacent to the central city, with a large area of urban waterfront parks, moderate intensity of construction and land development, and more reasonable planning and layout of blue and green spaces.With the benefits of urban construction and ecological services, the cooling

Spatial Distribution Patterns of LST
The heat island intensity is divided into five grades by mean standard deviation method.It uses the combination of average temperature and multiples of different standard deviations to define an urban heat island, which can accurately reflect the change in heat island intensity and avoid the difference in different phases [53].Figures 3 and 4 show that the portion of heat island patches in Jinan is 29.92%, mainly distributed in the central western and northeastern areas of the capital city, and the portion of cold island patches is 28.47%, distributed primarily in Daming Lake Park and Huashan Scenic Area.The hot island patches take up 28.68% of Qingdao, distributed in the south of Jimo District and the east of Chengyang District, while cold island patches take 28.46%, mainly distributed in the ecological, scenic areas such as Fufeng Mountain and Shimei Temple Park.The share of hot island patches in Heze is the highest, taking 31.19% of the city, distributed in the central urban area and the west of Mudan District, while the proportion of cold island patches is 29.61%, distributed in the northeast and southeast regions of the urban area.In terms of spatial distribution, the heat island areas are concentrated in the central urban areas where the proportion of impervious surface is high and the vegetation coverage is low, especially in the industrial and commercially packed areas with high population density, complex built environments, and scattered the green space patches.However, cold islands are mainly located in the suburbs and outer suburbs.The suburban area is adjacent to the central city, with a large area of urban waterfront parks, moderate intensity of construction and land development, and more reasonable planning and layout of blue and green spaces.With the benefits of urban construction and ecological services, the cooling effect of blue and green space is more noticeable.The outer suburbs are located at the edge of the built-up zone, with large forested scenic areas, extremely low population and building density, and a significant degree of cooling.The results show that high aggregation and connectivity of large blue-green patches have a more prominent cooling effect, while the highly aggregated impervious surfaces have a more prominent warming effect.
osphere 2023, 14, x FOR PEER REVIEW 9 of effect of blue and green space is more noticeable.The outer suburbs are located at the ed of the built-up zone, with large forested scenic areas, extremely low population and bui ing density, and a significant degree of cooling.The results show that high aggregati and connectivity of large blue-green patches have a more prominent cooling effect, wh the highly aggregated impervious surfaces have a more prominent warming effect.

Spatial Distribution Patterns of CID and CIS
The comprehensive assessment value of the CIS and the CID of the study area is vided into five grades by mean standard deviation method.As shown in Figure 5, t distribution of the CIS and the CID has evident spatial heterogeneity, and the high-val regions of CIS are mainly distributed in forest scenic areas at the edge of the urban bu up zone with extensive vegetation and sparse population.In contrast, the high-value C regions are primarily located in the core urban areas, especially in the dense industr and commercial zones and high-density residential areas.The low-value CIS regions a mainly located in areas with low green space coverage and high population density, wh the low-value CID regions are distributed in remote areas with low impervious surfa ratios and population densities.From the spatial distribution, it can be seen that the hig value CIS regions show the characteristics of multipoint distribution and independe dispersion, while the distribution of high-value areas of CID trends toward connectin gathering, and spreading around, with the high-demand and low-supply areas gradua overlapping.With the city's rapid development, the blue-green space is gradually occ pied and divided by urban buildings, the proportion of impervious surface keeps incre ing, and the supply capacity and demand level of the cold island become progressive unbalanced.sphere 2023, 14, x FOR PEER REVIEW 9 of effect of blue and green space is more noticeable.The outer suburbs are located at the ed of the built-up zone, with large forested scenic areas, extremely low population and buil ing density, and a significant degree of cooling.The results show that high aggregati and connectivity of large blue-green patches have a more prominent cooling effect, wh the highly aggregated impervious surfaces have a more prominent warming effect.

Spatial Distribution Patterns of CID and CIS
The comprehensive assessment value of the CIS and the CID of the study area is vided into five grades by mean standard deviation method.As shown in Figure 5, t distribution of the CIS and the CID has evident spatial heterogeneity, and the high-val regions of CIS are mainly distributed in forest scenic areas at the edge of the urban bu up zone with extensive vegetation and sparse population.In contrast, the high-value C regions are primarily located in the core urban areas, especially in the dense industr and commercial zones and high-density residential areas.The low-value CIS regions a mainly located in areas with low green space coverage and high population density, wh the low-value CID regions are distributed in remote areas with low impervious surfa ratios and population densities.From the spatial distribution, it can be seen that the hig value CIS regions show the characteristics of multipoint distribution and independe dispersion, while the distribution of high-value areas of CID trends toward connectin gathering, and spreading around, with the high-demand and low-supply areas gradua overlapping.With the city's rapid development, the blue-green space is gradually occ pied and divided by urban buildings, the proportion of impervious surface keeps increa ing, and the supply capacity and demand level of the cold island become progressive unbalanced.

Spatial Distribution Patterns of CID and CIS
The comprehensive assessment value of the CIS and the CID of the study area is divided into five grades by mean standard deviation method.As shown in Figure 5, the distribution of the CIS and the CID has evident spatial heterogeneity, and the high-value regions of CIS are mainly distributed in forest scenic areas at the edge of the urban built-up zone with extensive vegetation and sparse population.In contrast, the high-value CID regions are primarily located in the core urban areas, especially in the dense industrial and commercial zones and high-density residential areas.The low-value CIS regions are mainly located in areas with low green space coverage and high population density, while the low-value CID regions are distributed in remote areas with low impervious surface ratios and population densities.From the spatial distribution, it can be seen that the high-value CIS regions show the characteristics of multipoint distribution and independent dispersion, while the distribution of high-value areas of CID trends toward connecting, gathering, and spreading around, with the high-demand and low-supply areas gradually overlapping.With the city's rapid development, the blue-green space is gradually occupied and divided by urban buildings, the proportion of impervious surface keeps increasing, and the supply capacity and demand level of the cold island become progressively unbalanced.

The Correlation Coefficient and the Bivariate Moran's I between the CIS and CID
The statistical correlation coefficient of the CIS and the CID and bivariate glob ran's I results (Table 3) show that the correlation coefficient R is negative in all three indicating that the CIS and the CID are significantly negatively correlated and tha is an imbalance in matching the supply and demand for cooling in cold islands.high-value area of CID, because of the large population, the high density of bu alters wind turbulence and impedes ventilation and heat dissipation.Meanwhile, t pervious surface replaces natural surfaces, such as green spaces and water bodie their cooling effects.Therefore, CIS values tend to be low, and the bivariate global M I is negative.This finding indicates that the spatial distribution of the CIS and th showed an overall strong spatial dependency.That is to say, the increase in the local areas leads to a decrease in the CIS in the surrounding areas, and areas wit CIS values are adjacent to areas with low CID values, while areas with lower CIS t be close to areas with higher CID.This aggregation phenomenon indicates a sign spatial spillover effect between the CIS and the CID.   3) show that the correlation coefficient R is negative in all three cities, indicating that the CIS and the CID are significantly negatively correlated and that there is an imbalance in matching the supply and demand for cooling in cold islands.In the high-value area of CID, because of the large population, the high density of buildings alters wind turbulence and impedes ventilation and heat dissipation.Meanwhile, the impervious surface replaces natural surfaces, such as green spaces and water bodies, and their cooling effects.Therefore, CIS values tend to be low, and the bivariate global Moran's I is negative.This finding indicates that the spatial distribution of the CIS and the CID showed an overall strong spatial dependency.That is to say, the increase in the CID in local areas leads to a decrease in the CIS in the surrounding areas, and areas with high CIS values are adjacent to areas with low CID values, while areas with lower CIS tend to be close to areas with higher CID.This aggregation phenomenon indicates a significant spatial spillover effect between the CIS and the CID.

Matching Supply and Demand of CIS and CID
The degree of CIS and CID matching is characterized on the basis of the Z−score standardization and quadrant division method (Figure 6) noted in Section 3.5.There are four types of cooling effect matching of supply and demand: quadrant I (high supply−high demand), quadrant II (low supply-high demand), quadrant III (low supply-low demand), and quadrant IV (high supply-low demand).By calculating the proportions of the four matching states, as shown in Table 4, it can be seen that the matching of supply and demand for the cooling effect in the three cities is dominated by a state of low supply and high demand, strongest in Qingdao, followed by Jinan, then Heze.In all three cities, it accounts for more than 40% of their boundaries, indicating that most areas are facing the problem of insufficient cooling effect.The second type is high supply−low demand.More than 30% of Jinan's and the combined 20% of Qingdao's and Heze's territory are regions where the efficient utilization of the blue−green space cooling effect is not high.Finally, the matching states of high supply-high demand and low supply−low demand account for the smallest percentage of the three cities, only about 30%.The results show that the urban blue−green space cooling effect supply and demand capacities do not match, pointing to the problems of uneven distribution.Because of the unreasonable organization of blue−green spaces and the rapid expansion of urban construction, more areas may face the problem of an insufficient supply of cooling islands.Because of the great potential of blue and green space layout within the low supply-high demand areas, the government needs to strictly control the development intensity of urban construction, improve the layout of blue and green space, and increase the supply capacity of cooling effect.For areas with high supply and low demand, the government needs to adequately strengthen urban construction, reasonably organize blue−green space coverage, and improve the efficiency of cooling effect utilization.
Atmosphere 2023, 14, x FOR PEER REVIEW 11 of 18 the four matching states, as shown in Table 4, it can be seen that the matching of supply and demand for the cooling effect in the three cities is dominated by a state of low supply and high demand, strongest in Qingdao, followed by Jinan, then Heze.In all three cities, it accounts for more than 40% of their boundaries, indicating that most areas are facing the problem of insufficient cooling effect.The second type is high supply−low demand.More than 30% of Jinan's and the combined 20% of Qingdao's and Heze's territory are regions where the efficient utilization of the blue−green space cooling effect is not high.Finally, the matching states of high supply-high demand and low supply−low demand account for the smallest percentage of the three cities, only about 30%.The results show that the urban blue−green space cooling effect supply and demand capacities do not match, pointing to the problems of uneven distribution.Because of the unreasonable organization of blue−green spaces and the rapid expansion of urban construction, more areas may face the problem of an insufficient supply of cooling islands.Because of the great potential of blue and green space layout within the low supply-high demand areas, the government needs to strictly control the development intensity of urban construction, improve the layout of blue and green space, and increase the supply capacity of cooling effect.For areas with high supply and low demand, the government needs to adequately strengthen urban construction, reasonably organize blue−green space coverage, and improve the efficiency of cooling effect utilization.(1) The high-supply-low-demand areas are mainly distributed in scenic water areas and forest parks with sparse population and high green space coverage, such as Daming Lake Park and Qianfo Mountain Scenic Area in Jinan, Fufeng Mountain and Shimei Temple Ecological Park in Qingdao, and Huandi River Coast in Heze.The study [54] shows that large green space patches can enhance the infiltration and connection with the surrounding ecological patches (green spaces, water bodies, etc.) to a certain extent when the shape index increases, providing a good supply of cool island cooling capacity.Daming Lake Scenic Area is a large park green zone within the urban space, with a systematic and complete green ecological network connected inside and outside.The built environment of the scenic area is partly joint with the blue and green spaces where the water bodies and vegetation are coordinated and organized, which is more effective than the multiple superposition effect of the cooling capacity of a single water system and green space.Therefore, although located in the core  (4) The low-supply and low-demand zones are mainly distributed in the suburban adjacent to the edge of the central city, where the CIS is lower because of th island effect of the surrounding built-up areas.At the same time, the populat these areas, such as the edges of forest parks and scenic areas, is relatively low intensity of urban construction and development is moderate, and the CID is l On the basis of the coupling coordination model, the degree of coordinated dev ment between the CID and the CIS is calculated, as shown in Table 5.The average D in Jinan is 0.54, ranging from "extreme incoordination" to "good coordination".The age D value in Qingdao is 0.47, ranging from "extreme incoordination" to "good c nation".The average D value for Heze is 0.50, ranging from "severe incoordinatio "good coordination".Overall, the two primary development types of these cities are sitional development and failure and decline, with coordinated development taki more than one-third.Among the cities, Jinan has the highest proportion of coordi development type, followed by Heze, then Qingdao.Qingdao has the highest prop of transitional development and failure and decline types, followed by Heze, then It is worth noting that although the ecological livability index of Qingdao is highe that of Heze, its D-value ranks the lowest among the three cities, with a failure and d (1) The high-supply-low-demand areas are mainly distributed in scenic water areas and forest parks with sparse population and high green space coverage, such as Daming Lake Park and Qianfo Mountain Scenic Area in Jinan, Fufeng Mountain and Shimei Temple Ecological Park in Qingdao, and Huandi River Coast in Heze.The study [54] shows that large green space patches can enhance the infiltration and connection with the surrounding ecological patches (green spaces, water bodies, etc.) to a certain extent when the shape index increases, providing a good supply of cool island cooling capacity.Daming Lake Scenic Area is a large park green zone within the urban space, with a systematic and complete green ecological network connected inside and outside.The built environment of the scenic area is partly joint with the blue and green spaces where the water bodies and vegetation are coordinated and organized, which is more effective than the multiple superposition effect of the cooling capacity of a single water system and green space.Therefore, although located in the core urban area with a complex building population, Daming Lake Scenic Area still has a high cooling capacity of a cold island.(2) The low-supply-high-demand areas are mainly distributed in the central urban areas with a complex built environment, dense population, and low green space coverage, such as Tianqiao District and Huaiyin District in the west of Jinan and Shibei District near Jiaozhou Bay in Qingdao.Because of the high proportion of impervious surfaces, complex building structures, low albedo, and high imperviousness of human-made structures and surfaces, such areas easily absorb more solar radiation and reduce the amount of water used for evaporation and heat absorption, resulting in lower latent heat fluxes and increased latent heat fluxes.At the same time, these areas are located in commercial and industrial areas with dense population distribution and frequent socioeconomic activities, so the demand for cooling capacity of the cooling island is very high; the CID is much higher than the CIS, and the supply and demand matching of the cooling effect is extremely unbalanced.(3) The high-supply-high-demand areas are mainly distributed in the built-up areas adjacent to large parks and other green spaces, such as in Jinan High-tech District, Qingdao Shinnan District, and the southeastern part of Heze Mudan District.The building and population densities in those regions are relatively high, so the demand for the cooling effect is also high.These areas have a higher level of economic development and more financial resources to mitigate the heat island effect.Moreover, they are adjacent to large patches of green space, which can provide a certain amount of cold island supply because of the cooling "spillover effect" of the surrounding green area.Thus, the CIS value in those regions is higher.(4) The low-supply and low-demand zones are mainly distributed in the suburban areas adjacent to the edge of the central city, where the CIS is lower because of the heat island effect of the surrounding built-up areas.At the same time, the population in these areas, such as the edges of forest parks and scenic areas, is relatively low, the intensity of urban construction and development is moderate, and the CID is lower.

The D Value between the CIS and the CID
On the basis of the coupling coordination model, the degree of coordinated development between the CID and the CIS is calculated, as shown in Table 5.The average D value in Jinan is 0.54, ranging from "extreme incoordination" to "good coordination".The average D value in Qingdao is 0.47, ranging from "extreme incoordination" to "good coordination".The average D value for Heze is 0.50, ranging from "severe incoordination" to "good coordination".Overall, the two primary development types of these cities are transitional development and failure and decline, with coordinated development taking no more than one-third.Among the cities, Jinan has the highest proportion of coordinated development type, followed by Heze, then Qingdao.Qingdao has the highest proportion of transitional development and failure and decline types, followed by Heze, then Jinan.It is worth noting that although the ecological livability index of Qingdao is higher than that of Heze, its D-value ranks the lowest among the three cities, with a failure and decline type at over 20%.This may be due to the fact that Qingdao, as a coastal area, has a significantly higher level of urbanization development than inland areas, with noticeable differences in regional expansion, higher expansion compactness, and a lower degree of integrated and coordinated development.In addition, the supply of blue-green spatial elements in some regions has not been fully considered during the expansion, and the spatial layout of the artificially built infrastructure and blue-green environment is unreasonable, resulting in a low coupling between the supply and demand of cold islands.As seen in Figure 8, the lower D-value areas are mainly distributed in the low-coupling coordination areas where the positive interaction between the CID and the CIS is weak, such as the high-supplylow-demand regions in Huashan Scenic Area in the northern part of the central city of Jinan, Shimei Temple Park Scenic Area in the northern part of Licang District of Qingdao, the low-supply-high-demand areas in Tianqiao District and Shizhong District of Jinan, Chengyang District of Qingdao, and the eastern part of Mudan District of Heze.In these areas, the building density and blue-green space composition lack purposeful planning and coordination, and the supply capacity of the cold island effect is unbalanced.The high D-value areas are mainly concentrated in the high-supply and high-demand areas such as Jinan High-tech District, Qingdao Shibei District, and the western part of Heze Mudan District.The remarkable characteristics of those regions are relatively high building density and economic development level adjacent to the functional green park space with a relatively reasonable spatial layout.The supply capacity for the cold island effect is high, as is the level of demand for it.The supply and demand relationship between the CIS and the CID is balanced.Therefore, the ecological benefits provided by continuing to increase the blue-green infrastructure layout in the region are limited; thus, other ways to mitigate the heat island effect could be taken.Studies show that the uneven distribution of urban buildings, population, and blue-green space reduces the coupling and coordination of the CIS and the CID.This reminds us that in urban construction, it is necessary to thoughtfully lay out the blue-green space structure, especially in high-density built-up areas with low supply and high demand.Doing so may adequately reduce the intensity of urban development and construction, optimize the blue-green space network system, and form cold island cooling corridors and patches of different levels and types.This will ensure better matching of cooling supply and demand, resulting in their improved, coordinated development, and lead to better supply of the cool island effects of blue and green spaces.work system, and form cold island cooling corridors and patches of different levels an types.This will ensure better matching of cooling supply and demand, resulting in the improved, coordinated development, and lead to better supply of the cool island effec of blue and green spaces.

CIS and CID Comprehensive Evaluation Index System Establishment
The heat balance mechanisms within built-up areas are complex, and the built environment is often considered a significant heat source, while UGS can effectively cool the surrounding environment.However, the supply and demand relationship regarding the cooling effect of blue-green space is not yet clear, so it is necessary to explore how much cooling can be provided by UGS and how much cooling is needed to mitigate urban heat islands.We constructed the evaluation indicator system to quantify the cold island supply (CIS) and the cold island demand (CID) and revealed the law of spatial distribution and changes in the CIS and the CID.The CIS reflects the cold island supply capacity to mitigate the urban thermal environment.The CID demonstrates the level of cold island demand to mitigate the urban thermal environment.This quantitative evaluation is the essential precondition for identifying where the CIS and the CID need to be adjusted to achieve a regional equilibrium of supply and demand of the cooling effect of blue-green spaces.Therefore, from the supply and demand perspective, quantifying the combined level of public demand for cooling effects and the supply capacity of blue-green infrastructure can help identify key areas of imbalance between supply and demand for cooling effects, thus enabling the accurate allocation of blue-green infrastructure resources.

Spatial Distribution between CIS and CID
Exploring the spatial distribution patterns of the CIS and the CID allows for accurately identifying potential urban development and resource allocation problems [55].We have constructed a system of indicators for cold island supply capacity (CIS) and cold island demand levels (CID).On the basis of the results of the ecological livability assessment, three representative cities in Shandong were selected for testing.The results of the study show that there is significant spatial heterogeneity between the CIS and the CID.The high-value areas of CIS are mainly distributed in forest scenic areas at the edge of urban built-up areas with extensive vegetation and sparse population, with multipoint distribution and independent dispersion.The high-value CID regions are mainly located in dense industrial and commercial spaces and high-density residential areas, with a tendency to connect, cluster, and spread around.The reason for this is that along with the city's rapid development, the blue and green spaces are gradually occupied and divided by urban buildings, resulting in a gradual imbalance between the supply and demand capacities of the cold islands.This spatial distribution pattern indicates an uneven and contradictory distribution between the CIS and the CID.Therefore, clarifying the distribution characteristics and spatial patterns between the CIS and the CID can help to analyze the spatial configuration of matching supply and demand for cooling effects.

Spatial Relationship between CIS and CID
Exploring the spatial relationship between the CIS and the CID is important to optimizing the supply-demand matching relationship for the cooling effect in blue-green spaces [56].The correlation coefficient R and the bivariate Moran's I can reveal the spatial response law of the CIS and the CID.The results show that the CIS and the CID are significantly negatively correlated, with higher CID values tending to have lower CIS.There is an imbalance in the supply-demand match for cold island cooling and a significant spatial spillover effect.The coupled coordination model and the quadrant division method are used to analyze the supply and demand matching state of the cooling effect.The results show that the type of matching of the blue-green space cooling effect is a predominantly low supply and high demand, with most regions facing a spatial mismatch and quantitative imbalance of cooling effects.At the same time, the level of coupled and coordinated development of the cooling effect is not ideal, and the type of development is dominated by transitional development and failure to decline.Analyzing the spatial relationship between the CIS and the CID can help us develop targeted intervention strategies in urban planning.For example, in areas of high supply and low demand, we could make full use of cooling resources and optimize the urban spatial structure.In areas of low supply and high demand, development intensity should be strictly controlled, the proportion of blue-green spaces should be increased, and the design of blue-green spaces should be optimized.Within high-supply-high-demand areas, the continued increase in blue-green space provides limited ecological benefits, so other means should be used to mitigate the thermal environment.Therefore, exploring the spatial relationship between the CIS and the CID can help clarify the cooling effect of the blue-green space supply and demand relationship and provide a scientific basis for urban development and management.

Limitations and Prospects
This study focuses on the spatial relationship between the CIS and the CID and the changing law of matching supply and demand, which has certain theoretical and practical significance for urban development, construction, and mitigation of the heat island effect.However, the study is still in its infancy and has some limitations.First, the cooling effect supply and demand evaluation system is not comprehensive enough with respect to the selection of indicators and needs further improvement.There may also be a relationship between the CID and residents' age and type of employment that has not been included [57].For example, older people can be more sensitive to the heat island effect and need more cooling.The CIS may be related to the height and type of vegetation [58].Second, the mechanisms of flow change and future trends between the CIS and the CID have not yet been discussed, which imposes certain limitations on clarifying the supply and demand for blue-green space cooling effects.Third, research has focused on the assessment of cooling effects at the scale of built-up urban areas and the exploration of spatial relationships; however, the supply and demand for cooling effects of blue-green space at smaller or larger spatial scales have not been explored.Therefore, future research should focus on (1) establishing a more comprehensive evaluation index system for the CIS and the CID to accurately quantify the supply and demand for the cooling effect of blue-green space, (2) exploring the flow of trends between the CIS and the CID, and (3) conducting research on the matching of blue and green space cooling effects supply and demand at multiple spatial scales.

Figure 1 .
Figure 1.The study area of Shandong province, China.

Figure 1 .
Figure 1.The study area of Shandong province, China.

Figure 5 .
Figure 5. CIS and CID spatial distribution map: (a1) Jinan CIS, (b1) Qingdao CIS, (c1) Heze CIS, (a2) Jinan CID, (b2) Qingdao CID, and (c2) Heze CID.4.4.Spatial Relationship between the CIS and the CID 4.4.1.The Correlation Coefficient and the Bivariate Moran's I between the CIS and the CID The statistical correlation coefficient of the CIS and the CID and bivariate global Moran's I results (Table3) show that the correlation coefficient R is negative in all three cities, indicating that the CIS and the CID are significantly negatively correlated and that there is an imbalance in matching the supply and demand for cooling in cold islands.In the high-value area of CID, because of the large population, the high density of buildings alters wind turbulence and impedes ventilation and heat dissipation.Meanwhile, the impervious surface replaces natural surfaces, such as green spaces and water bodies, and their cooling effects.Therefore, CIS values tend to be low, and the bivariate global Moran's I is negative.This finding indicates that the spatial distribution of the CIS and the CID showed an overall strong spatial dependency.That is to say, the increase in the CID in local areas leads to a decrease in the CIS in the surrounding areas, and areas with high CIS values are adjacent to areas with low CID values, while areas with lower CIS tend to be close to areas with higher CID.This aggregation phenomenon indicates a significant spatial spillover effect between the CIS and the CID.

Figure 8 .
Figure 8. Distribution of the CCD between CIS and CID (a) Jinan (b) Qingdao (c) Heze.

Figure 8 .
Figure 8. Distribution of the D between CIS and CID (a) Jinan (b) Qingdao (c) Heze.

Table 1 .
Cold island supply capacity (CIS) index system.

Table 2 .
Cold island demand capacity (CID) index system.

Table 3 .
CIS and CID correlation coefficient and bivariate Moran's I statistics.

Table 3 .
CIS and CID correlation coefficient and bivariate Moran's I statistics.

Table 4 .
CIS and CID match type ratio statistics.

Table 4 .
CIS and CID match type ratio statistics.

Table 5 .
CIS and CID coupling coordination development ratio statistics.

Table 5 .
CIS and CID coupling coordination development ratio statistics.