Assessment of Urban Green Space Based on Bio-Energy Landscape Connectivity: A Case Study on Tongzhou District in Beijing, China

: Green infrastructure is one of the key components that provides critical ecosystems services in urban areas, such as regulating services (temperature regulation, noise reduction, air puriﬁcation), and cultural services (recreation, aesthetic beneﬁts), but due to rapid urbanization, many environmental impacts associated with the decline of green space have emerged and are rarely been evaluated integrally and promptly. The Chinese government is building a new city as the sub-center of the capital in Tongzhou District, Beijing, China. A series of policies have been implemented to increase the size of green urban areas. To support this land-use decision-making process and achieve a sustainable development strategy, accurate assessments of green space are required. In the current study, using land-use data and environmental parameters, we assessed the urban green space in the case study area. The bio-energy and its ﬂuxes, landscape connectivity, as well as related ecosystem services were estimated using a novel approach, the PANDORA model. These results show that (1) in the highly urbanized area, green space is decreasing in reaction to urbanization, and landscape fragmentation is ubiquitous; (2) the river ecology network is a critical part for ecosystem services and landscape connectivity; and (3) the alternative non-green patches to be changed to urban, urban patches which can improve landscape quality the most by being changed to green, and conservation priority patches for biodiversity purposes of urban green were explicitly identiﬁed. Conclusively, our results depict the spatial distribution, ﬂuxes, and evolution of bio-energy, as well as the conservation prioritization of green space. Our methods can be applied by urban planners and ecologists, which can help decision-makers achieve a sustainable development strategy in these rapidly urbanizing areas worldwide.


Introduction
Due to the increasing urbanization, as well as its direct and indirect impacts on biodiversity and ecosystem services, sustainable development strategy plays a decisive role in the urban areas around the world [1]. With the continued urban development, about 50% of the world's population and approximately 76% in developed countries are urban dwellers [2,3]. Besides, it is predicted that by the end of the 21st century, 90% of the world's population will live in cities [4,5]. Therefore, the ecological perspective of sustainable urban development is moving in the direction of a win-win strategy [6]. That is, we should simultaneously achieve ecosystem stability and human well-being in those places where most people live and can directly experience nature [7]. Some specific interventions are implemented in the urban ecosystem to achieve this win-win strategy, primarily through the optimized and reasonable assessment of green space [8,9].

Data Preparation
Multiple datasets were collected and compiled from diverse sources (Table 1). Forest location, forest types, forest conservation types, vegetation coverage, agricultural area location, and water body location were extracted from the Forest Inventory Database in 2015 for the Tongzhou District. The artificial surfaces land classes in 2015, such as buildings, roads, and railways, were obtained from a land-use map, which were interpreted from a series of Landsat Multispectral Scanner images and Landsat Thematic Mapper by Beijing Qianfan Shijing Technology Company (Beijing, China). The soil map was also extracted from the Forest Inventory Database. All administrative data, such as the boundaries of state, city, county (district), town (sub-district) and village, were provided by Beijing Municipal Bureau of Land and Resources Tongzhou Branch.  (Figure 1). In accordance with Corine Land Cover Classification System, land-use types were reclassified into 13 categories: continuous and dense urban fabric, roads networks, airports, permanently irrigated arable land, nurseries in permanently irrigated land, broad-leaved forest, coniferous forest, mixed forest, natural grasslands, moors and heathlands, bare rocks, sparsely vegetated areas, as well as rivers and streams ( Table 2).

Overall Framework
The main goal of the study was to evaluate the green space of Tongzhou District. An open-source plugin, PANDORA 3.0 model [35], was used to accomplish the assessment. This model ran in QGIS 2.16.3 software [36]. In this way, four indicators describing the landscape connectivity and ecosystem services values of green space were assessed.

1.
The bio-energy index (Mcal/year) describes the energy which an ecological system has to dissipate in the environment to maintain its level of metastability [37].

2.
To identify the bio-energy landscape unit (BELU) level of landscape connectivity, the bio-energy fluxes (Mcal/year) exchanged among adjacent BELUs was calculated. A BELU is an higher hierarchical level than the landscape mosaic, defined as a portion of landscape with variable homogeneity characteristics surrounded by recognizable and significant barriers for bio-energy fluxes [28,37].

3.
The dMtot index (%) indicates the connectivity importance of each patch. That is, a patch with a high dMtot value makes an enormous contribution of landscape connectivity to the overall ecosystem.
PANDORA is fundamentally a landscape evolution model based on a system of ordinary differential equations (ODEs) [18,[37][38][39][40]. An algebraic hierarchy solution was applied to simulate the final asymptotic bio-energy. Areas of high asymptotic bio-energy values indicated the most important regions for landscape connectivity [18]. On the contrary, those having low asymptotic bio-energy contributed less landscape connectivity. The Supplementary Materials provides a complete description and specific parameters for PANDORA. Also, see Pelorosso et al. [18,28] for more details of the model.
The theoretical framework of our study ( Figure 2) included four steps: (1) subdivision of landscape units; (2) bio-energy graph; (3) landscape evolution model; and (4) NUAs assessment at the patch scale. 2. To identify the bio-energy landscape unit (BELU) level of landscape connectivity, the bio-energy fluxes (Mcal/year) exchanged among adjacent BELUs was calculated. A BELU is an higher hierarchical level than the landscape mosaic, defined as a portion of landscape with variable homogeneity characteristics surrounded by recognizable and significant barriers for bio-energy fluxes [28,37]. 3. The dMtot index (%) indicates the connectivity importance of each patch. That is, a patch with a high dMtot value makes an enormous contribution of landscape connectivity to the overall ecosystem.
PANDORA is fundamentally a landscape evolution model based on a system of ordinary differential equations (ODEs) [18,[37][38][39][40]. An algebraic hierarchy solution was applied to simulate the final asymptotic bio-energy. Areas of high asymptotic bio-energy values indicated the most important regions for landscape connectivity [18]. On the contrary, those having low asymptotic bioenergy contributed less landscape connectivity. The Supplementary Materials provides a complete description and specific parameters for PANDORA. Also, see Pelorosso et al. [18,28] for more details of the model.
The theoretical framework of our study ( Figure 2) included four steps: (1) subdivision of landscape units; (2) bio-energy graph; (3) landscape evolution model; and (4) NUAs assessment at the patch scale.

Bio-Energy Graph
Subsequently, the bio-energy graph was illustrated in terms of bio-energy (see Equations (1) and (2), Supplementary Materials) and its fluxes index (see Equation (9), Supplementary Materials). The

Bio-Energy Graph
Subsequently, the bio-energy graph was illustrated in terms of bio-energy (see Equations (1) and (2), Supplementary Materials) and its fluxes index (see Equation (9), Supplementary Materials). The bio-energy graph consisted of nodes and arcs. Nodes were circles having dimensions proportional to the bio-energy of each BELU. The size of the circle was proportional to the bio-energy. Arcs were lines describing the fluxes of bio-energy between neighboring BELUs, whose thickness was proportional to the interaction of bio-energy. The bio-energy graph was available for graphically identifying the energy spatial distribution and fluxion of green space. Additionally, the final asymptotic values of each patch (see Equations (10) and (11), Supplementary Materials) were obtained to build an asymptotic bio-energy graph.

Alternative Scenarios Development
In order to detect the best land-use scale for urban planners, we conducted three alternative scenarios representing three expanding scales ( Figure 3). Scenario A (Figure 3a) is the sub-center area of Beijing proposed by the Urban Master Plan of Beijing (2016-2035) [31], and we adjusted the eastern boundary to the border of BELUs. Scenario B (Figure 3b) is the aggregation of the sub-center district and its adjacent towns. Scenario C (Figure 3c) refers to the entire administrative area of Tongzhou District. Then a comparison between different scenarios based on M as tot value (see Equations (8) and (12) in Supplementary Materials for further details) was conducted. M as tot value was the total asymptotic bio-energy of the overall system. Therefore, the evolution of M as tot index described the overall ecological quality of landscape connectivity. These evolution processes were compared among the three scenarios to regard the best land-use solution to be chosen. bio-energy graph consisted of nodes and arcs. Nodes were circles having dimensions proportional to the bio-energy of each BELU. The size of the circle was proportional to the bio-energy. Arcs were lines describing the fluxes of bio-energy between neighboring BELUs, whose thickness was proportional to the interaction of bio-energy. The bio-energy graph was available for graphically identifying the energy spatial distribution and fluxion of green space. Additionally, the final asymptotic values of each patch (see Equations (10) and (11), Supplementary Materials) were obtained to build an asymptotic bio-energy graph.

Alternative Scenarios Development
In order to detect the best land-use scale for urban planners, we conducted three alternative scenarios representing three expanding scales ( Figure 3). Scenario A (Figure 3a) is the sub-center area of Beijing proposed by the Urban Master Plan of Beijing (2016-2035) [31], and we adjusted the eastern boundary to the border of BELUs. Scenario B (Figure 3b) is the aggregation of the sub-center district and its adjacent towns. Scenario C (Figure 3c) refers to the entire administrative area of Tongzhou District. Then a comparison between different scenarios based on M as tot value (see Equations (8) and (12) in Supplementary Materials for further details) was conducted. M as tot value was the total asymptotic bio-energy of the overall system. Therefore, the evolution of M as tot index described the overall ecological quality of landscape connectivity. These evolution processes were compared among the three scenarios to regard the best land-use solution to be chosen.

NUAs Assessment
For the purpose of prioritizing the conservation and restoration area, an assessment of 835 NUA patches was implemented. The 835 NUA patches were alternative areas for afforestation by the local government. Beijing Tongzhou Forestry Bureau provided the location of the NUAs. We ranked the dMtot (see Equation (13), Supplementary Materials) and ESV_B values (see Equation (14), Supplementary Materials) of the selected NUA patches. Then these patches were reclassified into three types in terms of conservation prioritization. First, the areas where both dMtot and ESV_B were zero, were called Low Priority Zones. Second, the patches with a zero dMtot index but a nonzero ESV_B index were defined as Medium Priority Zones. Third, the other regions where neither dMtot nor ESV_B was zero were named as High Priority Zones.

Spatial Dynamics of Bio-Energy
The bio-energy graph ( Figure 4) shows the spatial dynamics of bio-energy obtained for the case study area. Nodes in Figure 4a display the normalized bio-energy of each BELU. Overall, the total bio-energy was 3.21 billion Mcal/year in 2015. Most of the high-level bio-energy BELUs (>10.0 Mcal/m 2 /year) were located in the northwest corner of the case study area, while the BELUs with lower bio-energy values were situated in the central areas (e.g., the sub-center area in Figure 3). Arcs in Figure 4a represent bio-energy exchange between neighboring BELUs. In general, 10.3% of the total bio-energy interacted between adjacent BELUs to maintain their landscape connectivity. The high-level bio-energy BELUs also represented high-levels of bio-energy flux. Meanwhile, the BELU 139 (marked by a red circle in Figure 4a) was the major part of the river ecology network, which associated with other 54 BELUs accounting for 42.9% of the total area. This result emphasized that BELUs of the river ecology network were more important for the landscape connectivity compared with other parts of the city ecosystem.

Detecting the Best Land-Use Scenarios
Land-use types and bio-energy evolution processes between three scenarios were compared to detect the best land-use scale. These three scenarios were adopted, representing three urban planning scales (Figure 3). Scenario A (Figure 3a) denoted a downtown scale with an area of 143.87 km 2 . Scenario B (Figure 3b) referred to an intermediate scale with an area of 647.99 km 2 . Scenario C (Figure 3c) stood for a landscape scale with an area of 906.5 km 2 . The expanding among Scenarios A to C expressed the urban sprawl process traditionally adopted by planning practice. Figure 5 shows the land-use variations in response to urbanization among three planning scales. The downtown Scenario A was the most urbanized area, which had a proportion of urban areas as high as 69.95% and the green The result of asymptotic bio-energy values of patches is demonstrated in Figure 4b. Colors of the patches mean asymptotic values, and the more essential patches for landscape connectivity exhibit deeper color [28]. Most high asymptotic bio-energy areas (>5.2 Mcal/m 2 /year) were also concentrated in the river ecology network. The representation highlighted that BELUs of the river ecology network were bio-energy centers closely connecting with the remaining BELUs. Comparatively, the other areas with low asymptotic bio-energy values contributed less landscape connectivity to the whole ecosystem.
In summary, we depicted the spatial distribution, fluxes, and evolution of bio-energy, as well as landscape connectivity. The bio-energy and asymptotic bio-energy graph (Figure 4) both identified that the river ecology network was a significant location for landscape connectivity, which deserves priority conservation. In contrast, landscape fragmentation was ubiquitous in the low asymptotic bio-energy areas in response to urban sprawl.

Detecting the Best Land-Use Scenarios
Land-use types and bio-energy evolution processes between three scenarios were compared to detect the best land-use scale. These three scenarios were adopted, representing three urban planning scales ( Figure 3 (Figure 3c) stood for a landscape scale with an area of 906.5 km 2 . The expanding among Scenarios A to C expressed the urban sprawl process traditionally adopted by planning practice. Figure 5 shows the land-use variations in response to urbanization among three planning scales. The downtown Scenario A was the most urbanized area, which had a proportion of urban areas as high as 69.95% and the green space only accounting for 18.71%. In contrast, the urbanization ratios of Scenarios B and C were lower than that of Scenario A, and the proportions of green space of Scenarios B and C were higher than that of Scenario A. These spatial urban expansion processes between Scenario A to C indicated that green space decreased in reaction to urbanization.   However, the bio-energy evolution in terms of M as tot index ( Figure 6) indicated a positive process trend for Scenario A, which had higher standardized M as tot values during the asymptotic process compared with Scenarios B and C. These results indicated that Scenario A had the best overall ecological quality of landscape connectivity. Comparatively, Scenario B had the worst evolution process with relatively low standardized M as tot values. Thus, though Scenario A was a high urbanization area, Scenario A had a better ecological quality of landscape connectivity than Scenarios B and C. Instead, Scenario B was identified to be a critical area for restoration. Some policies aimed at reasonably limiting urban sprawl should be conducted out of the sub-center area. Additionally, some conservation and restoration measures should be planned in the area of Scenario B to improve landscape connectivity.

Identification and Priority Ranking of Critical Areas from NUAs
To identify the critical areas for prioritizing conservation and restoration, 835 alternative afforestation patches were selected and evaluated. NUAs assessment based on dMtot and ESV_B indexes is illustrated in Figure 7. With dMtot ranking (Figure 7a Figure  7b) were scattered on the southeast side Beijing-Hangzhou Grand Canal. These results indicated that the river ecological network was the critical area for landscape connectivity and ecosystem services due to its bio-physical characteristics.
Then we reclassified the 835 NUAs in terms of dMtot and ESV_B indexes into three types ( Figure 8).
1. Low Priority Zones were the zero ESV_B areas. These zones consisted of artificial surface patches in NUAs. For example, the patches numbered 42, 68, 76, 323, 324, and 326 in Figure 8 are some buildings in a park. High Priority Zones surrounded most of the zones. Low Priority Zones should be reasonably utilized by inhabitants living around here. Current non-green patches to be alternatively changed to urban can be detected from these zones. 2. Medium Priority Zones with zero dMtot values but nonzero ESV_B values were caused by the isolation of urban fabric and major road networks. For instance, No. 358 and 311 patches in Figure 8 are cut off by the Sixth Ring Road thus lost the connectivity factor (dMtot index). Most In summary, we broadly detected the best land-use scale for planners among three urban sprawl scenarios. The comparison between different scenarios implied that rapid urbanization has caused green space decreasing. Furthermore, the downtown scale, Scenario A, was the best land-use scale. Some policies aimed at reasonable limiting urban sprawl and green space renewal projects should be conducted in the area out of Scenario A.

Identification and Priority Ranking of Critical Areas from NUAs
To identify the critical areas for prioritizing conservation and restoration, 835 alternative afforestation patches were selected and evaluated. NUAs assessment based on dMtot and ESV_B indexes is illustrated in Figure 7. With dMtot ranking (Figure 7a), the top ten patch numbers were 412, 27, 820, 176, 301, 331, 51, 270, 254, and 265. Almost all high dMtot level patches (deeper red in Figure 7a) concentratedly distributed on the banks of the Beijing-Hangzhou Grand Canal. While the patches numbered 412, 27,820,176,419,197,560,301,331, and 243 were sequentially ranked as the top ten ESV_B values (Figure 7b). The patches of high ecosystem services values (deeper red in Figure 7b) were scattered on the southeast side Beijing-Hangzhou Grand Canal. These results indicated that the river ecological network was the critical area for landscape connectivity and ecosystem services due to its bio-physical characteristics.
Then we reclassified the 835 NUAs in terms of dMtot and ESV_B indexes into three types (Figure 8).

1.
Low Priority Zones were the zero ESV_B areas. These zones consisted of artificial surface patches in NUAs. For example, the patches numbered 42, 68, 76, 323, 324, and 326 in Figure 8 are some buildings in a park. High Priority Zones surrounded most of the zones. Low Priority Zones should be reasonably utilized by inhabitants living around here. Current non-green patches to be alternatively changed to urban can be detected from these zones.

2.
Medium Priority Zones with zero dMtot values but nonzero ESV_B values were caused by the isolation of urban fabric and major road networks. For instance, No. 358 and 311 patches in Figure 8 are cut off by the Sixth Ring Road thus lost the connectivity factor (dMtot index). Most districts of the zones were isolated patches. To increase landscape connectivity within Medium Priority Zones, we should add some new elements like ecological corridors to the zones in current urban patches surrounding these zones, which can best improve the landscape quality.

3.
High Priority Zones were regions where green space was concentrated and contiguous. From a nature conservation and sustainability perspective, High Priority Zones were found to be the most important areas devoting to the overall city ecosystem. Some well-designed parks were included in the zones, such as Canal Forest Park in Figure 8. These zones were the priority protected districts for decision-makers. In the case of urbanization and related land-use change programs (e.g., road building, house construction) occurring in these zones, accurate environmental assessment should be implemented to avoid decreasing landscape connectivity.
Sustainability 2019, 11, x FOR PEER REVIEW 10 of 15 districts of the zones were isolated patches. To increase landscape connectivity within Medium Priority Zones, we should add some new elements like ecological corridors to the zones in current urban patches surrounding these zones, which can best improve the landscape quality. 3. High Priority Zones were regions where green space was concentrated and contiguous. From a nature conservation and sustainability perspective, High Priority Zones were found to be the most important areas devoting to the overall city ecosystem. Some well-designed parks were included in the zones, such as Canal Forest Park in Figure 8. These zones were the priority protected districts for decision-makers. In the case of urbanization and related land-use change programs (e.g., road building, house construction) occurring in these zones, accurate environmental assessment should be implemented to avoid decreasing landscape connectivity.

Landscape Fragmentation Caused by Urbanization
Landscape fragmentation, one of the consequences of urban expansion [41], is due to the patches of lost landscape connectivity through the spread of artificial surfaces [42]. Landscape connectivity can be measured not only by the characteristics of the landscape (structural connectivity), but also by the organism mobility (functional connectivity) [43]. Lots of studies have figured out that urbanization has caused the results of landscape fragmentation functionally and structurally [44][45][46]. For instance, landscape fragmentation resulting from urban expansion took place in Haidian District, Beijing, China [21]. In the current study, our results suggest that landscape fragmentation (low asymptotic bio-energy patches in Figure 4b) is occurring in our case study area as a result of rapid urbanization. Although green space covered 29.2% of Tongzhou District (Table 2 and Figure 5), most patches presented relatively low structural landscape connectivity due to isolation by artificial surfaces.

Spatial Reduction of Green Space in Response to Urbanization
Many studies have demonstrated that green space always decreases in response to urbanization [47][48][49]. Here, green space generally belongs to a natural or semi-natural open space with considerable vegetation coverage [50], including woodland (forest), grassland, and farmland [51]. For example, many Chinese cities, such as Dalian [49,52], Kunming [47], and Jinghong [23], have suffered temporal and spatial patterns of green space reduction in response to urbanization. In our study, the spatial relationship between the urban area and green space is broadly consistent with these studies. Figure 5 shows the spatial reduction of green space in response to urbanization among three planning scales. Among Scenarios A to C, the highly urbanized Scenario A generally had lower levels of urban greening than Scenarios B and C. Meanwhile, the relatively large bio-energy nodes in Figure 4a were generally distributed outside the downtown area. These results both suggest that urban development may reduce green space.

Landscape Fragmentation Caused by Urbanization
Landscape fragmentation, one of the consequences of urban expansion [41], is due to the patches of lost landscape connectivity through the spread of artificial surfaces [42]. Landscape connectivity can be measured not only by the characteristics of the landscape (structural connectivity), but also by the organism mobility (functional connectivity) [43]. Lots of studies have figured out that urbanization has caused the results of landscape fragmentation functionally and structurally [44][45][46]. For instance, landscape fragmentation resulting from urban expansion took place in Haidian District, Beijing, China [21]. In the current study, our results suggest that landscape fragmentation (low asymptotic bio-energy patches in Figure 4b) is occurring in our case study area as a result of rapid urbanization. Although green space covered 29.2% of Tongzhou District (Table 2 and Figure 5), most patches presented relatively low structural landscape connectivity due to isolation by artificial surfaces.

Spatial Reduction of Green Space in Response to Urbanization
Many studies have demonstrated that green space always decreases in response to urbanization [47][48][49]. Here, green space generally belongs to a natural or semi-natural open space with considerable vegetation coverage [50], including woodland (forest), grassland, and farmland [51]. For example, many Chinese cities, such as Dalian [49,52], Kunming [47], and Jinghong [23], have suffered temporal and spatial patterns of green space reduction in response to urbanization. In our study, the spatial relationship between the urban area and green space is broadly consistent with these studies. Figure 5 shows the spatial reduction of green space in response to urbanization among three planning scales. Among Scenarios A to C, the highly urbanized Scenario A generally had lower levels of urban greening than Scenarios B and C. Meanwhile, the relatively large bio-energy nodes in Figure 4a were generally distributed outside the downtown area. These results both suggest that urban development may reduce green space.

The Importance of Green Space Beside the River Network
Rivers in urban areas play an essential role in ecosystem services of many city ecosystems [18,28,53,54]. In this paper, the estimation of asymptotic values for bio-energy (Figure 4b) suggests BELUs of the river ecology network are also the most important ecological belt in the ecosystem. In addition, High Priority Zones in Figure 8 are concentrated on the sides of the Beijing-Hangzhou Grand Canal. Thus, both the asymptotic bio-energy map and the prioritization conservation map reveal that the river ecology network is a critical part for ecosystem services and bio-energy landscape connectivity in the case study area. However, the area of green space in the river ecology network represents only 7.95% of the green space among the overall ecosystem. Thus, these critical areas deserve careful environmental assessment in case of urban development or relevant land-use change programs [18]. Moreover, an increase of landscape connectivity of other green spaces to the river ecology network should be designed for the future development of urban forestry projects.

Contributions to the Literature
The current study applied and improved a new method, the PANDORA model, to evaluate urban green space based on bio-energy landscape connectivity. We believe that the viewpoint of landscape connectivity can efficiently promote studies and projects related to urban planning and biodiversity conservation at different spatiotemporal scales. Although other methods, such as FRAGSTATS (University of Massachusetts Amherst, Amherst, USA) [55] and Conefor Sensinode software (University of Lleida, Lleida, Spain) [56], alternatively provide arithmetic for calculation of landscape connectivity, PANDORA model solves the problem of data scarcity, using data that are usually available to land managers (e.g., remote sensing images, digital elevation models). Compared with other ecosystem services models (e.g., the InVEST [24], SPAN [25], and MIMES models [26]), PANDORA focuses on bio-energy interactions and landscape connectivity between patches. Urban planners could recognize whether the landscape fragmentation resulting from urbanization could take place in their cities. Ecologists can apply the model to ecosystem services related to landscape connectivity for biodiversity conservation purposes in urban contexts [18].

Conclusions
In this paper, we used an innovative approach, the PANDORA model, to evaluate the green space of Tongzhou District based on bio-energy landscape connectivity. We have achieved three main scientific goals (see the last sentence of the Section 1). Conclusions obtained from the three scientific goals show that: