Urbanizing the floodplain: global changes of imperviousness in flood-prone areas

Cities have historically developed close to rivers and coasts, increasing human exposure to flooding. That exposure is exacerbated by changes in climate and population, and by urban encroachment on floodplains. Although the mechanisms of how urbanization affects flooding are relatively well understood, there have been limited efforts to assess the magnitude of floodplain encroachment globally and how it has changed in both space and time. Highly resolved global datasets of both flood hazard and changes in urban area from 1985 to 2015 are now available, enabling the reconstruction of the history of floodplain encroachment at high spatial resolutions. Here we show that the urbanized area in floodplains that have an average probability of flooding of 1/100 years, has almost doubled since 1985. Further, the rate of urban expansion into these floodplains increased by a factor of 1.5 after the year 2000. We also find that urbanization rates were highest in the most hazardous areas of floodplains, with population growth in these urban floodplains suggesting an accompanying increase in population density. These results reveal the scope, trajectory and extent of global floodplain encroachment. With tangible implications for flood risk management, these data could be directly used with integrated models to assess adaptation pathways for urban flooding.


Introduction
More than half (55%) of the world's population lives in urban areas and that percentage is projected to increase to 68% by 2050 [1]. Cities are centers of wealth, infrastructure and economic activity, but also host impoverished and vulnerable communities. Any disruption to urban areas caused by flooding and other natural disasters therefore has the potential to impact human welfare [2]. The intensity of flooding is generally expected to increase globally with climate change [3][4][5]. At the same time, exposure to flood hazards is increased due to the urbanization of floodprone areas and the concomitant population growth. The effects of urbanization on floods are generally well understood [6][7][8], but it is also imperative to understand how urbanization has varied globally in relation to flood hazard.
Recent work shows that built-up areas on floodplains have increased [9], although most studies have been regional [10][11][12] with few global assessments of floodplain encroachment in both space and time [13,14]. A number of studies have examined the relationship between population dynamics, built-up areas and flooding at continental [15] and global scales [9,16,17]. Although these studies are tremendously valuable, their scope and methodological approaches presented a number of limitations in the context of assessing the urbanization of flood-prone areas. For example, the distinction between urban and non-urban areas on floodplains was addressed explicitly (e.g. [9,18]), implicitly (e.g. [16]) or not at all. Nonetheless, probably the most important factor for quantifying floodplain encroachment globally is the choice of datasets used to extract the pertinent information.
Global datasets of geolocated information on floodplains are not readily available with the exception of the recently developed GFPLAIN250m, which uses a morphometric approach to derive fluvial valley zones from elevation data [19]. Despite the global coverage of the GFPLAIN250m dataset though, it does have some limitations including its relatively coarse 250 m spatial resolution that could produce unreliable floodplain maps [20]. Moreover, the geomophologic approach used to derive GFPLAIN250m does not incorporate any of the physics governing water flow, while also not accounting for different levels of flood probability. An alternative approach could include mapping areas that became inundated from an adequately long-term dataset such as the Global Surface Water which classified 32 years of water surface dynamics from Landsat data [21], or the Global Flood Database that compiled MODerate resolution Imaging Spectroradiometer (MODIS) images from 2000 to 2018 and the Dartmouth Flood Observatory catalogue to derive flood event maps [17]. Nonetheless, the ephemeral nature of flood events, the inherent limitations of optical sensors (e.g. obscuring cloud cover) such as MODIS, and the difficulty in resolving flooding in urban areas imply that not all areas that could potentially inundate would be observed and information on flood probability would be difficult to extract.
In this paper, we attempt to overcome some of these limitations by combining maps of fluvial and pluvial flood inundation with maps of urban area changes to quantify floodplain encroachment globally as a function of flood probability. Although urbanization in flood-prone areas has increased overall, the magnitude, rates and associated hazard levels of floodplain encroachment have varied both in space and time and have not yet been systematically quantified globally.

Methods
We estimated floodplain encroachment by combining a dataset of annual urban maps with a model-derived floodplain delineation that assigns flood probability to ∼90 m pixels globally [22]. Both of these datasets offer better accuracy and higher spatial resolution compared with previously available datasets, improving our ability to resolve the localized impacts of flooding in urban areas on a global scale [23]. We then generated an annual time series of floodplain encroachment over the period 1985-2015 for different levels of flooding hazard by overlaying the urban area and floodplain maps (supplementary figure 1). In addition to estimating the area of floodplain encroachment, we determined the population potentially exposed to flooding of different probabilities given the urban sprawl on the floodplain.

Mapping changes in impervious area
The annual urban maps were obtained from the Global Annual Urban Dynamics (GAUD) dataset [24], which provides information on the year when urban expansion or reversal to increased vegetation (i.e. green recovery) occurred within the 1985-2015 period. GAUD is derived from Landsat imagery at 30 m resolution by fusing four global urbanization datasets (Global Human Settlement (GHSL) Layer, Global Urban Footprint, Global Urban Land, and the Global Artificial Impervious Area) and then applying the normalized urban areas composite index (NUACI) to identify years of urbanization or green recovery. The fusion algorithm kept the urban areas that were consistent among the other products, and used a locally adaptive classifier for the inconsistent pixels. The years of urban expansion and green recovery were identified with a temporal segmentation approach which fits a regression model to the NUACI time series and classifies the year with the largest residual as the year of change [25]. When validated against 200 000 samples, GAUD was shown to have higher accuracy compared to other available global datasets with lower commission and omission errors. Specifically, GAUD had an accuracy of detecting years of urbanization between 76% and 82% for the majority (90%) of global urban areas, while green recovery years were identified with an accuracy of 78% [24,26].

Calculating floodplain encroachment
We delineated floodplains globally using the Fathom-Global implementation of the LISFLOOD-FP flood inundation model [22], which is based on the local inertial formulation of the shallow water equations in 2D [27]. In order to account for different flood probability levels, the dataset contains ten discrete flood return periods including 5, 10, 20, 50, 75, 100, 200, 250, 500, and 1000 years for both fluvial and pluvial perils. River discharge for each of these simulations was estimated by a global regional flood frequency analysis [28], where the extreme flow behavior modeled in both gauged and ungauged catchments via transfer functions in the latter case. These flows were then routed in 1-D [29] through river channels (thresholded at 50 km 2 drainage area) delineated from the MERIT-Hydro hydrography dataset [30], with a conveyance constraint of a 2-year return period discharge that corresponds to bankfull flow conditions. Once bankfull capacity is exceeded, the model simulates water flow in 2D over the topography derived from the MERIT digital elevation model at 3 arc-second (∼90 m at the equator) spatial resolution. The pluvial model component generates nonriverine surface runoff via a rain-on-grid approach and simulates the inundation extent in catchments smaller than 50 km 2 , with extreme rainfall represented by Intensity-Duration-Frequency relationships [22,31]. Flood protection measures are estimated in the absence of a global inventory of such data. The FLOPROS database [32] is used to assign return period defense standards to river reaches, depending upon whether they are high-density urban, lowdensity urban, or suburban. Channels thus convey the 2-year discharge in rural areas, and the T-year discharge in other areas depending on the information in FLOPROS. Similarly, pluvial defense standards are assumed for different levels of urban density. For urban (high), urban (low), and suburban, assumed pluvial defenses correspond to the 2-, 5-, and 10-year in less developed regions and 5-, 10-, and 25-year in more developed regions.
Output of implementations of the model has been successfully and extensively validated in many regional case studies globally [31,33]. These studies show that the core numerical model can reproduce solutions to the full shallow water equations [27] and, when driven with gauged flows and high accuracy terrain data, can achieve accuracies of up to 90% (as determined using a robust Critical Success Index metric). When used with lower quality data, such as those available globally, these methods can still achieve Critical Success Index values of up to 70% at ∼90 m resolution, with errors in fractional inundated area decreasing below 10% when data are aggregated to 1 km 2 [22,23,31,34,35]. The flood maps generated from the model represent the maximum depth of pluvial and fluvial flooding in a 3 arc-second area during a flood event of a specific return period, and therefore were compressed into values corresponding to the highest frequency (i.e. lowest return period) flood when that pixel was inundated.
In order to quantify floodplain encroachment we overlaid the annual maps of urban area change onto the floodplain maps and then classified each pixel as encroached with an assigned flood probability if the pixel became urban during the study period  and was inundated for any of the flood simulations. Any areas where 'green recovery' occurred were reclassified as non-encroached floodplain starting at the year when recovery was detected according to GAUD. Flooding probability was defined from the return period of the simulation that was used to delineate the corresponding floodplain, but we also derived a flood hazard index which was calculated for a predefined area (0.1 • pixels). This aggregate hazard index was simply computed by multiplying the number of pixels (n f ) within the predefined area that were encroaching the floodplain for each return period, with the mean probability of occurrence ∑ where T f ∈ (5, 10, 20, 50, 75, 100, 200, 250, 500, 1000). The aggregate hazard index values along with the total encroachment area expansion for each 0.1 • pixel were binned into three classes each, with the bounds of each class derived from the 33th and 66th percentiles of the corresponding values, combining to form nine final classes.
There are a number of factors that can introduce uncertainty into the estimates of flood inundation (e.g. topography, boundary conditions) [36][37][38] and mapping of urban areas [39]. Since GAUD implicitly incorporates (via data fusion) multiple global urbanization datasets and the model-derived floodplain delineation is unique in terms of spatial coverage, resolution, and expicit accounting of flood exposure, we opted to use single datasets for mapping the floodplain and urbanization but report the results in a regional context (e.g. per country) to abate the effects of local-scale errors [40]. A comparison of the GAUD-derived floodplain encroachment with an independent reference dataset of urban imperviousness showed a relatively high agreement (supplementary figure 2).

Estimating population exposure to floods
Due to the variability of the available population datasets and the potential impact on estimates of population exposed to flood hazard [41], we used multiple datasets that satisfied the requirements of global coverage, spatial resolution less than 1 km, and multiple epochs. These datasets included: the 9 arc-second (∼250 m) GHSL dataset that represents population during 1990, 2000, and 2015 [42]; the Global Rural Urban Mapping Project that has 30 arc-second (∼1 km) resolution and represents 1990 and 2000 [43]; the Worldpop dataset that provides annual data from 2000 to 2020 on global population at a 3 arc-second spatial resolution [44]; the Gridded Population of the World dataset that represents population at 30 arc-seconds resolution every 5 years since 2000 [45]; and Landscan TM that provides annual population data since 2000 at 30 arc-seconds resolution globally [46]. Due to the scale discrepancy between the derived floodplain encroachment data (3 arc-seconds) and all the population datasets (9 and 30 arc-seconds) with the exception of Worldpop, we calculated the fractional areas of floodplain encroachment per flood probability as the number of 3 arc-second floodplain pixels for each return period T (masked to retain only pixels that were urbanized) divided by the number of non-null pixels within each 9 or 30 arc-second pixel. The population size for a given flood probability was then derived by multiplying the fractional area by the number of people in each 9 or 30 arc-second pixel. In the case of Worldpop, we just aligned the population and floodplain encroachment grids as their spatial resolutions matched resulting in a direct calculation of the population exposed to flooding.
We also calculated the growth rates for both urbanization and population for the entire and subsets of the study period as where x is the urban area or population for the ending and beginning years, t and b respectively. When assessing population changes the periods we used were 1990-2000, 2000-2015 and 1990-2015, whereas urbanization rates were calculated either for each year or the 1985-2000, 2000-2015 and 1985-2015 periods. We performed the same calculations for nonurban areas, that is areas that were not classified as urban during the specific period in the GAUD dataset. On average, the urban area with a 1/100 annual probability of flooding increased by 2.5% per year, and this rate of encroachment was relatively constant for the Americas, Europe and Africa (supplementary figure 3). Asia had an average 5.6% encroachment rate with a peak of 8.9% during the 30-year study period. These relatively high rates were also evident when assessing the proportion of urbanization that occurred after 1985 ( figure 1(b)). Almost all encroachment in the 100-year floodplain occurred after 1985 in India, Philippines, Thailand, Vietnam, Indonesia and parts of China as well as parts of the Middle East, sub-Saharan Africa, Spain, Mexico, Brazil, and United States. On the other hand, increased levels of urbanization in the 100-year floodplain for some regions like central Europe, western United States, Japan and southeastern China appears to have occurred prior to 1985 with little additional encroachment in the past three decades.

100-year floodplain encroachment
As well as quantifying floodplain encroachment as a fraction of the total urbanized area, we can also examine the associated population occupying those areas. As of 2015, 457.1 ± 134.7 million people in urban areas are exposed to floods with a 1% annual exceedance probability (i.e. 1/100-year floodplain). The large uncertainty in that estimate confirms results from previous studies that highlighted discrepancies among population datasets when estimating flood exposure [41,47]. Most of that population lives in Asia, however the proportion of global population in the floodplain is higher in Asia as well. Relative to their total population, the urban population residing in the 1/100-year floodplain varies substantially: 6.8 ± 2.6% for Asia, 8.0 ± 0.7 and 6.5 ± 2.2% for North and South America, 5.6 ± 1.6% for Europe, 3.6 ± 1.5% for Africa, and 3.8 ± 0.2% for Oceania. The 2015 population in both urban and non-urban areas that is potentially exposed to a 1/100-year flood is 1.49 ± 0.03 billion, which is consistent with recent studies [48].

Change in global floodplain urban areas
The proportion of urban areas exposed to flooding through time is shown in figure 2, which further breaks down the proportions by the highest frequency at which a given area is inundated. The urbanized area in the 1/100-year floodplain appears to have practically doubled between 1985 and 2015 (from 8.2% to 16.2%), which is also the case for all flood frequencies.
In addition to the increase of floodplain encroachment globally, the rate of annual urban expansion has increased since 2000 for all flood probabilities (p ⩽ 0.01). In fact, it appears that floodplain urban areas have expanded about 1.5 times faster on average during 2000-2015 compared to the 1985-2000 period. Floodplains corresponding to return periods less than 100-years had the highest rates of urbanization globally throughout the study period with average rates between 2.3 and 2.4% per year. Almost half of the 2015 urban area in the 1/100-year floodplain resides in the 1/20-year floodplain (i.e. on land with a higher flood frequency by a factor of 5). The urban expansion rate was slightly slower for floodplains with smaller annual flooding probabilities at 2.2% and 2.1% per year.
Encroachment of relatively high flood probability areas appears to have intensified since 1985 in many regions in Asia, Europe, and North America      (1990-2015, 1990-2000, and 2000-2015) in relation to flood annual probability (expressed in terms of the return period T).

Floodplain urbanization and population exposure
Finally, we assess global floodplain encroachment in conjunction with total population for different levels of flood hazard. The number of people in urban areas on floodplains (regardless of hazard level included here) increased from 0.45 (±0.21) billion in 1990 to 0.87 (±0.23) billion in 2015, while the corresponding numbers were 1.18 (±0.15) and 1.48 (±0. 19) respectively for non-urban areas. Despite the variability between the different datasets, there is consensus that the non-urban proportion of the population exposed to flooding is larger than the urban proportion ( figure 4(a)). About 3/4 of the non-urban population on the floodplain are exposed to flooding probabilities of 1/100 years, while the same is true for about half the urban population. These proportions are consistent among the different datasets, while also remaining relatively constant for the three epochs (1990, 2000, and 2015) available in the population datasets used. On the other hand, we found a significant difference between the distribution of the overall population increase on the floodplain between urban and non-urban areas. Specifically, 33.2 ± 0.8% of the population increase globally between 1990 and 2015 occurred on floodplains, and most of that increase was found in urban areas (20.7 ± 6.7% of the total population increase).
Total population on the floodplain has increased at a slower rate when comparing the 1990-2000 and 2000-2015 periods, with the highest rates associated with elevated flooding probabilities (return periods of 100-years or less). The slower population growth after 2000 was true for both urban and non-urban areas on floodplains, but the discrepancy between the rates was significantly more pronounced for non-urban areas ( figure 4(b)). More importantly, population growth was not only higher in urban areas compared to nonurban areas on floodplains but is also increased for higher flooding probabilities. The highest population growth rate in urban floodplain areas was 4.2 ± 0.6% for the 1/20-year return period and the lowest was 3.2 ± 0.4% for the 1/1000-year probability. Despite the variability in the datasets, we consistently found that the population size has increased in urban areas on the floodplain. It is difficult to attribute the accelerated urbanization rates on the floodplain as driven by population or related to commercial and industrial purposes [24], but these results do agree with recent work that found that built-up areas and population on the floodplain have increased in most of the world's countries [9,18].

Discussion and conclusions
Urbanization is driven by a number of factors including economics (e.g. job opportunities), industrialization, population pressure, political conditions, natural disasters in rural areas [49] and is expected to triple its 2000 extent by 2030 [50]. The link between urbanization and potential increase in flood risk stems from hydrologic alteration from impervious surface areas [51,52], effects on local climate [53], or changes in exposure (e.g. population distribution changes) [54]. Recent studies found that urbanization can significantly exacerbate flood risk; for example, during Hurricane Harvey the probability of flooding was elevated by ∼21 times due to urbanization [55]. Therefore, our results that showed the doubling of the urban area within the 1/100-year floodplain and the accelerated rate of urbanization since 2000 have significant implications despite some potential limitations.
Although the simulation that we used to delineate the floodplain globally relied on a hydrodynamic model, it only approximately accounted for flood protection measures suggesting that our results correspond to a potential 'worst case' (pluvial and fluvial) flood exposure. The presence of flood defenses such as dams and levees would reduce the flood risk to the exposed population [56], and there have been recent efforts to improve their representation in high-resolution, large-area models [57]. On the other hand, there have been a number of studies positing that the installation of structural protection measures changes the perception of risk at an institutional level leading to the further urbanization and population increase (or resettling after a flood event) in floodprone areas [16]. In the context of this 'safe development paradox' , the actual efficacy of flood defenses is further complicated when considering issues such as effects of levees [58], damages and maintenance requirements, or the socio-economic conditions of the population suffering flood losses or benefiting from protection measures [59,60]. Furthermore, even areas that are protected by defenses build to a specific standard (e.g. 100 years) can become exposed in the near future due to increases of flood magnitudes for the same return period [23,61].
In sum, we conclude that (i) human settlement on floodplains is increasing, with urbanization rates accelerating; (ii) the total urban area in floodplains has nearly doubled since 1985; (iii) population size has grown on the floodplain with the highest rates in areas with 1/20 annual probability of flooding. We argue that our analysis can act as the baseline for assessing the relationship between urbanization and flood hazard, as the resulting dataset is physically-based at sufficiently high spatial resolution with global coverage. Nonetheless, a number of open questions remain and can be addressed by expanding on this study. The hydraulic model simulations are dependent on the static DEM elevations, and therefore only represent a snapshot of flooding (circa 2000). As increases in impervious area as well as other changes in the landscape can affect flooding characteristics, allowing for dynamic floodplain delineation would be necessary to assess the feedbacks between urbanization and flood hazard globally and account for non-stationarity. Moreover, predicting future flood exposure due to the compounding effects of urbanization and population growth with expected changes in extreme rainfall return periods would similarly require dynamically coupled modeling, Finally, an integrated framework would be needed to better understand the socio-economic drivers of floodplain encroachment and inform the different aspects of national/local flood management.

Code availability
The code for the LISFLOOD-FP hydrodynamic model is available at https://doi.org/10.5281/zenodo. 4073011. Processing and visualization code is available upon request from K Andreadis.