Ecological redline policy may significantly alter urban expansion and affect surface runoff in the Beijing-Tianjin-Hebei megaregion of China

Urban expansion leads to surface changes that disrupt hydrological processes and increases flooding risks in cities. This increase may be severe in urban megaregions where clusters of cites have agglomerated. The China Ecological Redline Policy (ERP) is a national policy that protects priority areas with high-value ecosystem services. However, it is not clear how the ERP alters megaregion expansion and what this means for surface runoff across entire regions. By integrating specified models, we developed future urban expansion scenarios for 2030 with and without the ERP in the Chinese Beijing-Tianjin-Hebei (BTH) megaregion. The annual surface runoff volume under the ERP scenario decreased by 78 million m3 compared to the non-ERP involved scenario, but the ERP effectiveness at surface runoff regulation was different between the ecological redline areas (ERAs) and the non-ERAs. This suggested that multi-solutions should be incorporated into megaregions, such as regional ERPs and local, nature-based solutions, which could efficiently reduce the risk of urban flooding across whole regions.


Introduction
The development-environment conflict has become a major challenge across the globe, particularly in areas where rapid urbanization has increased the tension between urbanization, environmental protection, and economic growth (Bai et al 2012, Du et al 2018, Wing et al 2018. During the urban growth process, large areas of land are converted to impervious surfaces, which are known to affect surface runoff, and subsequent flood control and management in many urban areas (Trinh and Chui 2013, Hoekstra et al 2018, Zhang and Chui 2019. Increased runoff coefficients and shorter concentration times lead to increased flood risk (Nirupama and Simonovic 2007, Kjeldsen 2010, Guneralp et al 2015. Many cities have suffered urban flooding. For example, an urban flooding event that occurred in Beijing in 2012 caused 70 deaths and economic losses of 16 billion RMB . Furthermore, city agglomeration is an important economic development strategy in China. The government has promoted the emergence of large-scale urban agglomeration development modes. Many cities have now agglomerated and closely interact with each other. This has resulted in significant changes to land use patterns and surface runoff across a whole region rather than a single city (Day et al 2014, Sahana et al 2018. Therefore, it is very important to assess how the expansion of urban agglomeration affects surface runoff, and to explore possible measures and coordinated development in cities that could potentially improve flood mitigation. Despite a growing need to find policy solutions that minimize the negative impacts of urban expansion on surface runoff, most practical solutions remain limited to local scales. At a local scale, the main actions have focused on urban storm water management, such as the introduction of low-impact development plans in the United States, the sustainable urban drainage system in the United Kingdom, and the sponge city proposed by China (Hamel et al 2013, Mcgrane 2015, Satyavati and Shirishkumar 2018. These efforts are effective solutions for regulating urban surface runoff that utilize the ecological functions (e.g. infiltration, storage, and interception) of specific land uses/covers. These are also known as nature-based measures and include green roofs, storm water parks, and permeable pavements, etc (Elliott and Trowsdale 2007, Imran et al 2013, Zolch et al 2017, Akther et al 2018, Du et al 2020. However, few studies have focused on mitigation solutions at larger scales, such as urban megaregions. This is due to the gathered development patterns of urban megaregions. These development patterns have led to more complicated land use/cover dynamics and consequent surface runoff than the ones found in individual cities. Some studies have shown that optimizing the spatial layout of megaregion expansion should reduce the risk of urban flooding (Wiering et al 2017).
The Ecological Redline Policy (ERP) was produced by the Chinese national government and its aim was to resolve the conflicts between regional conservation and development (Jia et al 2018). This policy has established priority conservation zones through the use of ecological redline areas (ERAs) at a regional scale, and these redline areas are based on ecosystem services assessments (Bai et al 2018). Ecosystem services are potentially important to human well-being, and are deemed an effective approach for identifying ERAs (Zhang et al 2020). Currently, the ERP is being promoted at the national scale (Jia et al 2018), and several provinces or cities are currently delimitating their ecological redline areas (Bai et al 2016). Some studies have assessed the effects of the ERP on urban expansion (Jia et al 2018, Bai et al 2018. For example, Bai et al (2018) found that the ERP in Shanghai led to a decrease in constructed area and an increase in forest compared to the baseline scenario. However, there has been no attempt to determine the most efficient ways of integrating the ERP into megaregion expansion and investigate the subsequent effects on surface runoff.
In this study, we examined the impact of the ERP on urban expansion and its potential effect on surface runoff across the Beijing-Tianjin-Hebei (BTH) megaregion in China. This large urban region is one of the most important megaregions in the eastern part of north China and is facing a serious developmentenvironment conflict. The aims were to investigate (1) ways of integrating the ERP into megaregion expansion; (2) what effects the ERP would have on megaregion expansion; and (3) how the alterations to megaregion expansion by the ERP would further affect surface runoff in the megaregion. First, we simulated the future land use/cover patterns of the BTH megaregion under different scenarios. These were the regular trend scenario (RT), which did not consider the impact of the ERP, and the ecological protection scenario trend (EPT), which did include the ERP. Second, we simulated the annual surface runoff under the two different scenarios and analyzed the impacts of different megaregion expansion scenarios on the surface runoff process at regional to subregional scales. Finally, we briefly discuss the policy and scientific implications of the results.

Land use/cover (LULC) scenarios
We created three LULC scenarios in order to analyze the differences due to LULC change with and without the ERP and the further effects on surface runoff in the BTH region. These scenarios were Baseline 1 (BL) is a baseline scenario, which is the actual LULC map in 2015.
Scenario 1 (S1) is the regular trend (RT) scenario for 2030 without considering the impact of the ERP on LULC allocation.
Scenario 2 (S2) is an ecological protection trend (EPT) scenario for 2030, which considers the effect of the ERP.

Land use/cover classification
We used the Landsat-TM images from June-October, 2015 to generate the baseline land use/cover map (BL), with a spatial resolution of 30 m × 30 m. There are six land use/cover types in the BTH region, which are wetland, barren land, forest land, impervious land, grassland, and cultivated land. A backdating approach and an object-based method were integrated to identify the land use/cover types in 2015. We used the land use/cover classification of 2010 as the basic reference map to extract the changes in the land use/cover between 2010 and 2015, and then the change area was further classified. We used the stratified random sampling scheme to select locations within the cities and compared the image classifications to reference land use data created from the visual interpretation of high spatial resolution SPOT images (2.5 m) in eCognition Developer (www.ecognition.com/). The overall accuracy, Kappa statistic, and average producer and user accuracies exceeded 93%, 0.8, 80%, and 94%, respectively.

Evaluation of ecosystem services
A national ecosystem services assessment has been implemented in China (Ouyang et al 2016) and the assessment methods used in a previous study were applied to the BTH megaregion . Based on the importance of ecosystem services to this megaregion and the availability of spatial data, we selected six ecosystem services and then used the national ecosystem services assessment methods to map ecological redline areas in the BTH region for 2015. The six services were soil retention, water retention, carbon sequestration, sand storm prevention, flood mitigation, and habitat provision for biodiversity. The data sources and data reliabilities are shown as table 1. (1) Soil retention: this refers to soil reserved by ecosystems. The universal soil loss equation (USLE) was used to calculate soil retention, and was calculated by the InVEST model. The USLE can be expressed as follows (Yang et al 2016): where SC represents the soil retention capacity (t/ha), R is the rainfall erosivity factor, K is the soil erodibility factor, LS is the topographic factor, and C is the vegetation cover factor.
Among them, R was calculated by the daily rainfall erosivity model and the input data was the daily precipitation data, which was acquired from the China National Meteorological Information Center; K was closely related to soil characteristics and was estimated using the erosion/productivity impact calculator (EPIC); the soil data was obtained from the China Soil Science Database; LS reflected the impact of slope length and steepness on soil retention and was calculated based on an empirical formula (Wei et al 2002); and C was obtained from relevant foreign and domestic literature (Wei et al 2002, Rao et al 2014. (2) Water retention: This refers to water retained in ecosystems. Water retention was estimated using the water balance equation (Gong et al 2017): where TQ is total water retention (m 3 ); P is precipitation (mm); R is storm runoff (mm); ET is evapotranspiration (mm); and A is the area of a specific ecosystem (m 2 ).
Major data inputs for the water retention calculation include rainfall data, runoff data, and evapotranspiration data for 2015. The evapotranspiration data were obtained from the Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Science. The storm runoff was calculated by multiplying precipitation by the runoff coefficient, which was estimated from previous studies (Gong et al 2017). All the input data were standardized and we used the spatial analyst tool in ArcGIS to calculate the total water retention.
(3) Carbon storage: this refers to the carbon in terrestrial ecosystems. After referring to the national ecosystem services assessment (Ouyang et al 2016), we calculated the biomass carbon storage of three ecosystem types in the BTH region. These were forest, grassland, and wetland ecosystems. The calculation formula is as follows (Lewis et al 2009): where BCS is the biomass carbon storage; B i is the biomass density of ecosystemi(t/km 2 ), which was obtained from the Institute of Remote Sensing and Digital Earth, Chinese Academy of Science; CC i is the carbon content in the ecosystem biomass, which is 0.5 for forest and wetland, and 0.45 for grassland (Fang et al 2010); and AR i is the area of each grid cell.
(4) Sand storm prevention: This refers to sand reserved by ecosystems. We used the revised wind erosion equation (RWEQ) to estimate the sand storm prevention service, which was a function of the weather factor (WF), soil erodibility factor (EF), soil crust factor (SCF), surface roughness (K'), and vegetation cover parameter (C) (Borrelli et al 2017). The sand storm prevention volume was equal to the difference between the soil loss (S L ) under the bare land condition and the vegetation cover condition The RWEQ model can be expressed as follows: where Q max is the maximum transport capacity (kg/m); S is the critical field length (m); Z is the distance from the upwind edge of the field (m); and S L is the soil loss (kg/m 2 ). Among them, WF was calculated by dividing the total wind value for each period by 500 and multiplying by the number of days in the period (Ouyang et al 2016); the weather data were obtained from the China National Meteorological Information Center; EF reflected the relationship between wind erosion and soil physical and chemical properties; and SCF reflected the resistance of soil aggregates and crusts to windblown sand and was calculated using the SCF equation (Borrelli et al 2017). C has a significant effect on sand storm prevention and the vegetation cover data was acquired from the Chinese Academy of Sciences. The fraction of the ground surface covered with non-erodible plant material, and the silhouettes from standing plant residues and growing crop canopies were used in the RWEQ model to estimate the effect of vegetation cover (C) on sand storm prevention (Bilbro and Fryear 1994). The original RWEQ model was for the field scale, whereas our study was at a regional scale. Therefore, we used the roughness caused by the topography to replace the soil ridge roughness (K'), and was calculated using the Smith-Carson equation (Ouyang et al 2016).
(5) Flood mitigation: this refers to the flood mitigating volume of wetlands (e.g. lakes, reservoirs, and swamps). The national ecosystem services assessment (Ouyang et al 2016) was used to evaluate the flood mitigation capacity of various wetland types in the BTH megaregion. The area of wetlands was derived from the LULC data, and the storage capacity of wetlands was obtained from the Records for Chinese Lakes and The China Water Statistical Yearbook.
For lakes, available storage capacity was used as an indicator to measure flood mitigation capacity, which was a function of lake area (Ouyang et al 2016). The calculation formula is as follows: where C l is the available storage capacity (10 4 m 3 ) and A is the lake area (km 2 ). The coefficient was obtained from Ouyang et al (2016).
For reservoirs, the flood control storage capacity (C r ) was used as the flood mitigating indicator (10 4 m 3 ), which was calculated based on its total storage capacity ((10 4 m 3 )), and the coefficient was obtained from Ouyang et al (2016): For swamps, the surface stagnation of water (C s ) was used as the flood mitigation indicator (10 4 m 3 ). C s was calculated by multiplying the swamp area (A (km 2 )) by the average maximum depth in a flood period (D (cm)), and 1 m was used as the depth (Zhao et al 2003): (6) Provision of habitat for biodiversity: this refers to the total habitat area covered by indicator species. China performed a national biodiversity assessment in 2010 (Ouyang et al 2016). Due to the lack of species category and suitable habitat data, we adopted the corresponding evaluation results for the BTH region, which were based on the national assessment result (Ouyang et al 2016).
The national biodiversity assessment process was as follows: firstly, a total of 2820 species were selected as indicators to estimate biodiversity, including endemic, endangered, and nationally protected species. Then different protection targets (60%, 40%, or 30%) were assigned to species with different levels of rarity so that site selections could be carried out. Secondly, a habitat mapping process was developed using a simplified conceptual model (Xu et al 2009), and MARXAN (Klein et al 2010), a site selection software program, used an irreplaceability index to measure the biodiversity conservation value for each analytical unit. In our study, we used the 'high' and 'very high' areas of the national biodiversity assessment to delineate ecological redlines. The simplified conceptual model is as follows: where pH i indicates whether polygon i is a potential habitat or not; C i is the historical distributions of the indicator species in county i; I is the overlap of suitable elevation, slope, and aspect for the species; and H is suitable habitat for the species (Xu et al 2009).

Delineation of the ecological protection redline areas
We created an integrated index based on the six ecosystem services estimated above. This index, together with the information about ecosystem services and the relative abilities of different areas to supply them were used to map the ERAs.
The delineation of ERAs was undertaken using ArcGIS and the processes were as follows (Ouyang et al 2016): (1) we calculated the supply of each ecosystem service in a grid cell; (2) to quantify the relative importance of each grid cell for a single service, all grid cells were sorted by the supply of a specific ecosystem service in descending order, calculated the cumulative proportion of the specific ecosystem service across all grid cells, and then classified each grid into one of four levels of importance (very high, high, medium and normal). For example, we assigned 'very high' to the grids with cumulative proportion between 0% and 50%, 'high' to the grids with cumulative proportion between 50% and 75%, 'medium' to the grids with cumulative proportion between 75% and 90%, and 'normal' to the grids with cumulative proportion between 90% and 100%; (3) By overlay analysis, the relative importance of each service was synthesized into a composite index of ecosystem service importance using the maximum values method. The composite index value equaled the highest importance value of any ecosystem service in each grid cell. Finally (4), the 'very high' and 'high' grid cells were defined as the ERAs.

Land use/cover scenarios
We used the CLUE-S model to simulate the different land use scenarios for 2030 (S1 (RT) and S2 (ERT)) based on the ERAs. The CLUE-S model is a spatially explicit model that can integrate the analysis of land use/cover change with the driving forces. The model mainly had two main parts, a non-spatial demand module and a spatially explicit allocation module (Verburg et al 1999(Verburg et al , 2002. For the non-spatial module, the areal demands for different land use types on a yearly basis of the study area need to specify. Then, for the spatially explicit allocation module, the areal demands for different land use types will be allocated at specific location according to driving forces. The main advantages of this model were that it could simulate different scenarios and different land use/cover types at the same time by integrating spatial and nonspatial driving factors (Liu et al 2017a).
To simulate the LULC allocations in the BTH megaregion under the two different scenarios for 2030, in the RT scenario, we determined the yearly areal demands of different land use types based on the annual average LULC change from 2000 to 2015 by using the linear trend extrapolation .
In the EPT scenario, we determined the yearly areal demands of impervious land as similar to that in the RT scenario, while the areal demand of the wetlands, the forest, cultivated land a was according to the goals of the Wetland Conservation Plan (2015-2030) of Hebei, and the Master Plan (2006-2020) of Beijing, Tianjin and Hebei. Then the CLUE-S model allocated the areal land use demands for the time interval of one year. In the RT scenario, there was no spatial restrictions, and in the EPT scenario, the ecosystem redline areas were set as a spatial restriction. The logistic regression model was used to identify the potential factors driving the spatial allocation of land use/cover (Islam et al 2018). We selected 14 driving forces based on previous research (Liu et al 2017b). These 14 driving forces were divided into three classes, which were the geographical, social, and economic factors (table 3). To simplify the analysis, all the variables were transformed into the raster format (Verburg et al 2002). In addition, the goodness of fit of the logistic regression model was evaluated using the curve of receiver operating characteristic (ROC) method, which can evaluate the predicted probabilities by comparing them with the observed values over the whole domain (Phung et al 2019). The area below the ROC curve represents the ROC value, whose prediction accuracy improves as the value gets closer to 1 (Lamichhane and Shakya 2019). Finally, we used the CLUE-S model to generate future land use scenarios. The Kappa coefficient was used to assess the effectiveness of the CLUE-S model. It evaluates simulated probabilities by comparing them with ideal probabilities. In this study, the LULC simulation result for 2015 was compared with the interpreted LULC map to calibrate the land use type specific conversion and assess the accuracy of the CLUE-S model.

Simulation of average annual surface runoff
We estimated the average annual surface runoff under different land use scenarios using the long-term hydrological impact assessment (L-THIA) model, which was developed by Purdue University in the United States.
The L-THIA model is a distributed hydrological model, and can be used to simulate surface runoff on a large spatial-temporal scale (Lim et al 2006). The L-THIA model is based on the SCS-CN (Soil Conservation Service-Curve Number) method (Shadeed andAlmasri 2010, Li et al 2018), which simulates the runoff volume for each grid cell from the CN value and daily rainfall. The CN value is a dimensionless parameter that comprehensively reflects land use type, soil group, and antecedent moisture conditions (Bhaduri et al 2000).
The calculation formula for the L-THIA model is as follows: where, Q is runoff (mm); P is rainfall (mm); I a is initial abstraction (mm); and S is potential maximum retention after runoff begins. The data required for the model included land use data, soil data, hydrological data, and meteorological data. The soil data for the BTH megaregion were obtained from the China Soil Database and classified into four categories (A, B, C, and D) according to the classification standard of the United States Department of Agriculture, as required for L-THIA model. Soil permeability decreased from A to D. The daily runoff data from 1950 to 2015 were used to operate the L-THIA model and were acquired from The China National Meteorological Information Center. The L-THIA model was validated using observed runoff data from Luanxian reservoir, Wangkuai reservoir, and Xidayang reservoir, and the relative error was ± 10%, which indicated that the model could accurately simulate the surface runoff in the BTH region (Ju et al 2020) (figure 1). The CN values for the BTH region were obtained from Ju et al (2020) and are shown in table 2.
The surface runoff volume is closely related to the size of accumulated surface area. Therefore, we used the normalized average annual runoff depth (NAARD) to characterize the surface runoff generation ability by eliminating the impact of city size on surface runoff. The calculation formula is as follows (Chen et al 2017): where NAARD city is the normalized average annual runoff depth for a city (mm); RV city is the average annual runoff volume for a city (m 3 ); and A city is the administrative area of a city (m 2 ).

The ecological redline areas in the BTH region
In general, the ERAs covered about 106 173.36 km 2 in total and accounted for 49.53% of the BTH region (figure 2). The spatial distribution of the ERAs in the BTH region was uneven, and the northern and western mountains, the eastern Bohai Bay coastline, and inland lakes were the primary ecological reserves. The ecosystem services within the ERAs also had different spatial distributions (figure 2). Within the ERAs, the soil retention conservation areas were concentrated in the Yanshan and Taihangshan mountains in the northern and western parts of the BTH region. They covered 103 245.10 km 2 . Water retention zones within the ecological redline areas covered 72 277.81 km 2 . They were mainly distributed along the eastern coastlines, in inland lakes, and the boundary between mountains and plains. Carbon sequestration areas within the ERAs were mainly distributed in the Yanshan Mountains of the northern BTH region, and their area covered 69 405.00 km 2 in total. Sand storm prevention and conservation areas covered 59 637.20 km 2 , and only occurred in the northwestern forest of the BTH region. The flood mitigation area distribution was similar to water retention and the area within the redline area covered 81 510.50 km 2 . The biodiversity areas covered 31 845.79 km 2 , and mainly occurred in the northern, central, and southern parts of the BTH region.

Urban expansion of the BTH megaregion under the different scenarios
The logistic regression results are shown in table 3. The ROC value for each LULC type was greater than 0.7, which indicated that the regression model was suitable. Comparing the LULC simulation result with the interpreted LULC map for 2015 allowed us to obtain a Kappa coefficient of 0.79, which showed that the CLUE-S model performed satisfactorily.
The total BTH megaregion expansion patterns for 2030 under the different scenarios were similar, but the change ranges for the different LULC types varied ( figure 3, table 4). For example, the increases in impervious land across the BTH megaregion under the two 2030 scenarios were similar when compared to the baseline scenario in 2015, but their spatial occurrence was distinct. The impervious land in the BTH megaregion increased by 26.34% under the RT scenario relative to 2015 and increased by 25.02% under the EPT scenario. A comparison of the RT and EPT scenarios showed that there was an obvious difference in urban expansion inside and outside the ERAs. Inside the ERAs, the increase in impervious land was 1890 km 2 under the RT scenario for 2030 relative to 2015, but there was no increase under the EPT scenario. Outside the ERAs, the increase in impervious land was 411.04 km 2 under the RT scenario for 2030 relative to 2015, whereas the increase was 5707.08 km 2 under the EPT scenario. The increase in impervious land outside the ERAs under the EPT scenario for 2030 was about 1589.04 km 2 more than under the RT scenario relative to the baseline in 2015. For example, the increase in impervious land across Baoding and Handan under the EPT scenario was 94.32 and 104.04 km 2 more than under the RT scenario, respectively.
In addition, the cultivated land decreased under both 2030 scenarios, but the cultivated land under the EPT scenario decreased slightly less than that under the RT scenario. In contrast, under the EPT scenario, there was substantial increase in wetland relative to the RT scenario. For example, the total wetland area for 2030 under the EPT scenario increased by 9.20% compared to 2015, but decreased by 8.36% under the RT scenario. The forest area under the EPT scenario increased 221.04 km 2 more than under the RT scenario. The most obvious difference in land use change between the RT and EPT scenarios occurred in the northern part of the BTH megaregion, especially in Beijing, Chengde, and Zhangjiakou. In contrast, the land use changes in the other regions within the BTH megaregion were similar under the two future scenarios.  At the sub-region level, the impervious land in most cities across the BTH megaregion had expanded by 2030 relative to 2015 under both 2030 scenarios, but there were certain differences between the scenarios (table 5). Generally, there were two major differences between the two 2030 scenarios. In one case, urban expansion under the EPT scenario was less than under the RT scenario. There were three cities following this type of expansion, which were Beijing, Tianjin, and Chengde. The largest decrease in impervious land had occurred in the mega-city of Beijing by 2030 under the EPT scenario compared to the RT scenario, and the reduced area was 697.32 km 2 . There were ten cities where urban expansion under the EPT scenario for 2030 was faster than under the RT scenario. Among these cities, the largest increase in impervious land by 2030 under the EPT scenario occurred in Zhangjiakou, at 156.24 km 2 .

Impacts of BTH megaregion expansion on surface runoff under the different scenarios
Generally, on regional scale, the total average annual surface runoff for the entire BTH megaregion by 2030 under both the EPT and RT scenarios increased relative to that in 2015 (table 6). The total surface runoff for 2030 under the RT scenario increased by 5.87% (8.46 hundred million m 3 ) relative to 2015, and increased by 5.33% (7.68 hundred million m 3 ) under the EPT scenario. The megaregion expansion under the EPT scenario reduced potential surface runoff by 78 million m 3 across the entire BTH megaregion compared to the RT scenario. In addition, the contributions of the different land use/cover types to surface runoff across the BTH megaregion were different. What was noteworthy was that the surface runoff generated from impervious land made the largest contribution to the total runoff among all land use/cover types in 2030, whereas it was farmland that contributed the most to total runoff in 2015. The contribution to surface runoff made by impervious land was 32.83% in 2015, and will be 38.56% and 38.68% in 2030 under the RT and EPT scenarios, respectively.
There was an obvious difference in the surface runoff change between inside and outside the ERAs from 2015 to 2030 under the RT and EPT scenarios (table 7). The total surface runoff for 2030 under the EPT scenario decreased by 4.45% (2.57 hundred million m 3 ) relative to the RT scenario inside the ERAs, whereas it increased by 1.88% (1.78 hundred million m 3 ) outside the ERAs. Furthermore, in the two future scenarios, forest land made the largest contribution to the surface runoff inside the ERAs because it accounted for 50% of the total surface runoff. Impervious land and cultivated land had the two highest surface runoff contribution rates outside the ERAs, and the total contribution made by the two land use types exceed 90% under both future scenarios.
At the city scale, the normalized average annual runoff depth (NAARD) for most cities increased under both the RT and EPT scenarios by 2030 relative to that in 2015, but there were also distinct among cities for the two scenarios ( figure 4, table 8). Compared to the baseline scenario in 2015, Tianjin city showed the largest increase (17.71 mm) in NAARD, followed by Beijing (10.82 mm) under the RT scenario. In contrast, under the EPT scenario, Tianjin city also showed the largest increase in NAARD at 15.96 mm, but Beijing was ranked fourth (6.04 mm). Thus, the impacts of the ERP on the different cities in terms of city-level surface runoff were different to the RT scenario. There were three cities where the ERP could reduce the NAARD relative to the RT   '1' to '4' represents 'very high' , 'high' , 'medium' , and 'normal' , respectively. scenario. For example, the NAARD for Beijing by 2030 under the EPT scenario decreased by 4.79 mm (-5.52%) relative to that under the RT scenario. In contrast, the urban expansion under the EPT scenario might increase the NAARD of ten cities compared to the RT scenario. Among these cities, Handan showed the largest NAARD increase. The EPT results for Handan showed that there would be an additional 1.71 mm (1.75%) of NAARD relative to the RT scenario.

Discussion
Under the current development trends, megaregions in China will continue to expand, and this will inevitably intensify the conflicts between development and the environment (Fang et al 2018). This is a considerable challenge that decision makers must confront. Our results showed that under the RT scenario for 2030, surface runoff in the BTH would increase by 8.46 hundred million m 3 in total, compared to the  baseline in 2015, which is nearly 1/3 of total domestic water consumption in Beijing. This may reduce the already limited water resources and increase the flood risk that the BTH megaregion is currently facing . Therefore, researchers and policy makers in China need to minimize the negative impact of megaregion expansion on surface runoff. Our study attempted to explore whether the national ecological redline policy has an effect on surface runoff across a large urban megaregion in China,  and through scenario analysis, we found megaregion expansion under the EPT scenario can result in a lower surface runoff increase relative to the RT scenario. This finding revealed that altering urban expansion under the ERP scenario can protect the specific priority areas that have important ecological services, and reduce overall surface runoff generation by the megaregion. Although the surface runoff reduction (0.78 hundred million m 3 ) under the EPT scenario relative to RT scenario was limited, this study shows the potential value of the policy-making alternatives when attempting to promote the development of large-scale megaregions. Our findings also confirmed the results of previous studies that investigated these effects at individual city scales. For example, Bai et al (2012) reported that the inclusion of the ERP in urban land use planning in metropolitan Shanghai, China, increased terrestrial habitat protection by 174% and maintained flood mitigation. Our results also showed the different runoff changes inside and outside the ERAs under the EPT and RT scenarios. They showed that the impacts of the ERP on different sub-regions of the BTH megaregion were distinct. Although the surface runoff generated inside the ERAs under the EPT scenario was much less than that under the RT scenario, the greater surface runoff generation outside the ERAs under the EPT scenario may offset some of the surface runoff reduction within the ERAs relative to the RT scenario. This is because the inclusion of the ERP into the megaregion expansion plans prohibited the growth of impervious land within the ERAs, and thus, promoted an increase in impervious land outside the areas in order to maintain the same expansion rate for the entire region under the EPT scenario as would occur under the RT scenario. This leads to more potential surface runoff generation outside the ERAs under the EPT scenario than under the RT scenario. Therefore, these results emphasize that the implementation of the ERP alone into the spatial land use planning for megaregion expansion may be insufficient to reduce the surface runoff in all sub-regions within the megaregion. In the case of the BTH megaregion, cities, such as Handan and Baoding, which are not included in the ecological redline areas, need to have more urban expansion to compensate for the urban growth constraints in ERAs. Therefore, sub-regional solutions need to be created for these cities in order to regulate urban surface runoff. There are already some effective city-level solutions that mitigate surface runoff, such as low-impact development measures . Therefore, if joint regional-level ERP and sub-regional solutions are integrated into the spatial planning policies for megaregions, then surface runoff would be efficiently minimized in both regions and sub-regions.
Our study also had some limitations. For example, there is no standardized approach for identifying ERAs due to the diversity and complexities of ecosystem services. Furthermore, it was difficult to acquire high-quality data, and parameterize and interpret models etc. These limitations could lead to difficulties when attempting to identify potential ERAs. Therefore, there is a need to standardize procedures for these issues so that the identification of ERAs can become more effective. Additionally, the quality of the driving factors data used for the urban expansion simulation may have affected modeling accuracy. Furthermore, during model calibration (L-THIA model), the lack of available data might have affected the modeling reliability. We also did not consider the impacts of future climate change on surface runoff, which will play an important role in runoff generation. In further studies, alternative models and higher-quality data could be used to simulate urban expansion patterns and hydrological processes.

Conclusions
This study proposes a new way to examine the effectiveness of including the ERP in large-scale urban megaregions by analyzing how the ERP can alter urban expansion and effect the alteration on surface runoff. Through scenario analysis, the inclusion of the ERP into BTH megaregion expansion plans could reduce the annual surface runoff in the entire region while maintaining a similar growth rate for urban expansion envisaged under the RT scenario. We showed that the ERP effectively promoted urban expansion and regulated surface runoff in megaregions. These findings will significantly affect sustainable development, and urban flood control and management. However, the impacts of the ERP on surface runoff were distinct within and beyond the ERAs. The increase in surface runoff outside the ERAs may offset some of the surface runoff reduction when the EPT scenario is compared to the RT scenario. This implies that the ERP alone may not be able to mitigate surface runoff in all cities within megaregions. Integrating regional ERPs with beyond-ERA solutions can help minimize the effects of urban expansion on surface runoff at the sub-regional scale Our study has the potential to be scaled up to find modes of sustainable growth in urban megaregions throughout and beyond China.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.