Variation in the hydrological response within the Quebrada Seca watershed in Costa Rica resulting from an increase of urban land cover

ABSTRACT Urbanization is a global phenomenon which has provoked severe disruptions in hydrological cycles, resulting in flooding problems. While detailed studies exist for the world’s temperate zones, they are few for tropical zones where most of future urbanization may occur and where flooding is already a problem. A tropical watershed in Costa Rica was used to analyze the urban development and the associated hydrological response between 1945 and 2019, based on remotely sensed data and a numerical model. Using a detailed spatial-temporal approach, we found that the watershed’s overall urbanization over the timespan (+64%-points urban-areas) had led to major hydrological challenges (+80% runoff-volume, +220% peak-flow-rate and maximum-specific-discharge, and −25 min time-to-peak). These challenges were then placed in the context of historically reported flood events, providing a basis for spatially-differentiated flood mitigation actions and for guiding future urbanization. The study also provides valuable insights for other tropical regions with the same situation.


Introduction
More than half of the world's population already lives in urban areas. Forecasts predict that this share will rise to two-thirds by 2050 (United Nations 2018). The spatial distribution of megacities has also changed significantly: while the world's largest cities were previously located in richer countries in higher latitudes (e.g. London and New York), most urban growth is now taking place in developing countries in subtemperate and tropical zones (Day, Gunn, and Robert Burger 2021). Looking at environmental impacts, urbanization results inter alia in significant land cover changes, in turn severely changing hydrological responses. The latter aspect has been amplified by the construction of extensive stormwater drainage systems to convey rainwater out of urban areas.
Looking at the effects of urbanization on flood occurrence, there is consensus that increased impervious land cover reduces infiltration and thus increases runoff volumes and the magnitude of any floods (Villarini and Slater 2017), estimated effects on flood volumes range from zero (Dudley et al. 2001;Wright et al. 2013) to increases of significant magnitudes (Prosdocimi, Kjeldsen, and Miller 2015). Apparently, it is difficult to quantify the specific effects of changes in impervious land cover on flooding because many rivers are also responding to concurrent changes in climate, irrigation practices, land cover, water withdrawals and transfers, and flow regulation through dams (Allaire, Vogel, and Kroll 2015;Hejazi and Markus 2009;Viglione et al. 2016). Besides an increase in flood volumes, urbanization effects are also often related to flashier storm hydrographs, hence higher peak flow rates and lower times-topeak (Walsh et al. 2005).
While the impact of urbanization on flooding has largely been studied in countries with temperate climates, fewer studies deal with tropical countries, with the exception of Singapore (Chaosakul, Koottatep, and Irvine 2013). Most tropical countries are developing countries with high degrees of urbanization (United Nations 2018), much of which is unplanned and devoid of flood remediation measures (Dos-Santos et al. 2021). Moreover, there is often a lack of official information on the urbanization process, on the related infrastructure development, and on flood-related effects (Vaux et al. 2020). One way to overcome this lack of information is to use remotely sensed data to determine land cover change over time due to urbanization (Khan, Chapa, and Hack 2020). While aerial photographs (orthophotos) have been available since the early 1900s, they are often costly and not regularly taken, thus providing a low temporal resolution. However, multispectral satellite imagery with continuous high temporal resolution has been freely available for all places around the globe since the 1970s.
To assess the hydrological impacts of land cover changes resulting from urbanization, hydrological models are commonly used (Osuide 2021). One urban hydrological model used throughout the world for analyzing, planning, and designing stormwater runoff and drainage systems is the Storm Water Management Model (SWMM) of the United States Environmental Protection Agency (EPA). The model enables all hydrological processes to be spatially modelled by dividing a study area into a collection of smaller, homogeneous drainage areas (Rossman and Huber 2015). The degree of urban land cover can be represented in each of these areas by fractions of pervious and impervious sub-areas, while overland flow can be routed between sub-areas, between drainage areas, or between drainage system entry points. This makes SWMM suitable for the hydrological modelling of flood-related impacts due to urbanization (Rossman and Huber 2015).
Latin America has experienced rapid and unplanned process (Barros 2004), with populations in many countries concentrated in metropolitan areas subject to frequent urban flooding (Vaux et al. 2020). In this study, the spatial-temporal urbanization of the Quebrada Seca watershed in the Greater Metropolitan Area (GAM) in the Central Valley of Costa Rica is analyzed using remotely sensed data for a land cover classification between 1945 and 2019 for a total of 11 uniformly distributed years. The tropical watershed was chosen as our case study due to repeated flooding and a ruling of the Costa Rican constitutional court not to compound the flooding problems by further urbanization (Oreamuno and Villalobos 2015). The watershed is used to study urbanization in the GAM and the resulting change in the hydrological response. The urban hydrology software PCSWMM, which combines SWMM with a geographic information system (GIS) interface, was used with a spatial discretization into 17 elementary drainage areas, enabling a detailed spatial analysis of flood-and runoffgenerating processes over time.
The main objectives of this study are 1) to analyze long-term land cover change as a result of urbanization; and 2) to analyze flood generation due to land cover changes using a hydrological model. This enables the identification of the flood-producing effect of different degrees of urbanization in relation to antecedent land cover conditions and to put them in the context of historically reported flood events.
The specific novelty of this study is the combination of 1) a long-term land cover change analysis based on remotely sensed data and 2) the hydrological modelling of a rapidly urbanizing tropical watershed. Both the land cover change analysis and the hydrological modelling are characterized by a high spatial and temporal resolution. The detailed spatial identification of how runoff is generated and evolves temporally provides a meaningful basis for determining the location of effective source control measures. Hence, the results of this study support decision-making on the spatial distribution of flood prevention measures.
The study was conducted in the context of the SEE-URBAN-WATER (www.see-urban-water.uni-hannover.de) research project on the planning, design and implementation of naturebased solutions to reduce urban flooding while providing socio-ecological benefits.

Study area
Part of the San José Great Metropolitan Area (GAM), the 23km 2 Quebrada Seca watershed is located just northwest of the capital in the central valley of Costa Rica. Its main watercourses are the Quebrada Seca River, the Burío River and Aries Creek (Masís-Campos and Vargas-Picado 2014). Primarily of volcanic origin, the watershed covers parts of six municipalities, San Rafael de Heredia, Barva, Heredia, Flores, Belén, and Alajuela ( Figure 1). Its mean slope is 20%, with an elevation ranging from 869 to 1626 m.a.s.l. (Chen et al. 2021). Under natural conditions, the volcanic rock has developed into highly permeable soils (Masís-Campos and Vargas-Picado 2014) well-suited for recharging aquifers and controlling surface runoff. However, 66% of the total watershed area is now urbanized, directly impacting hydrological processes through reducing infiltration and increasing and accelerating surface runoff.
The climatic region in which the Quebrada Seca watershed is located has average annual precipitation of around 2000 mm (Solano and Villalobos 2012;Masís-Campos and Vargas-Picado 2014) and an average annual temperature of 24.8°C (Chaves et al. 2014). Located on the Pacific side of the country, the watershed has a well-defined bi-modal rainy season (Chen et al. 2021) featuring maximums in May-June and in August-September-October, separated by the midsummer drought in July (Maldonado et al. 2013).

Land cover classification
An analysis of the temporal variation of land cover in the Quebrada Seca watershed was carried out through classifying satellite images from Landsat sensors, taken between 1975 and 2019. For the two years, 1945 and 1960, where no satellite images were available, aerial images in panchromatic bands were obtained from National Geographic Institute of Costa Rica. For those years, mosaics were built from the aerial images using the software ERDAS Imagine 2015, with land cover classification performed manually.
For the period from 1975 to 2019, Landsat images were obtained from different sensors: Landsat 2 for 1975 to 1979, Landsat 5 for 1985 to 2001, Landsat 7 for 2007 to 2013 and Landsat 8 for 2019. The data was downloaded from the United States Geological Survey (USGS) Earth Explorer. Satellite images were considered in five-year intervals to avoid missing important changes in land cover. Images were selected for the months between December and April, corresponding to the dry season, a period where there is a better chance of gaining cloudless images (clouds can substantially impact classification accuracy). Table 1 provides detailed information on the aerial and satellite images used in this study.
A supervised classification based on the maximum likelihood methodology was used to determine land cover. Classification was performed with the extension 'Spatial Analyst' of ESRI ArcGIS 10.5. Land cover classes were defined based on field surveys, expert knowledge and previous work performed in the watershed, where the main land cover classes in different years were identified (Masís-Campos and Vargas-Picado 2014; Chaves et al. 2014;Oreamuno and Villalobos 2015). Four different classes were defined: high vegetation, low vegetation, pastures and urban cover. Each of these classes features different dominant land cover and uses, as detailed in Table 2. The classification was designed to reflect the study's hydrological approach, taking account of the effects of the historical changes in land cover on the watershed's hydrological response.
Classes used for the land cover classification and related, sub-summarized dominant land covers and uses

Rainfall runoff modelling
To evaluate the watershed's hydrological response due to ongoing urbanization, land cover distributions for the analyzed years were used as input for the hydrological model, while all other model parameters remained unchanged. Hydrological model was divided into 17 elementary drainage areas, in order to carry out a spatial analysis of the hydrological response to urbanization at different spatial scales: watershed, subwatershed, and elementary drainage area level. However, the hydrological modelling results could not be compared with actual streamflow measurements, due to the lack of hydrometric and channel geometry data for the specific event modelled.
Three sub-watersheds were defined using control points. The first control point was placed at the confluence of the Burío and Quebrada Seca rivers creating the Burío River subwatershed and the Quebrada Seca sub-watershed. This enabled an analysis of the headwaters in the upper part. The second control point was positioned at the intersection of the river with the highway on Route 1, a critical and flood-prone point, enabling an analysis of the hydrological response at this point  within the watershed, i.e. the sub-watershed at Highway 1. The third control point was located at the outlet of the whole watershed in the canton of Belén. Data from the Santa Lucía rain gauge located near its center ( Figure 1) was used to characterize precipitation events in the watershed. Annual series of maximum daily rainfall were acquired for the period 1982 to 2020 from the National Meteorological Institute database. An extreme event analysis was conducted to estimate expected precipitation for different return periods. Carried out with the 'extRemes' package of the R software (Gilleland and Katz 2016), the analysis was based on the Generalized Extreme Value Distribution method. To model the typical behavior of storms often occurring in the watershed, the characteristic temporal distribution for extreme events ( Figure A1 in Appendix A) as registered by the Santa Lucía rain gauge was used. The selected precipitation event has a volume of 130.9 mm and corresponds to an event of ten years of return period.
To simulate rainfall runoff, the software PCSWMM Version 7.3.3095 (Rossman and Huber 2015) was employed, a dynamic model combining the US EPA Storm Water Management Model (SWMM) with a geographic information system (GIS) interface. The hydraulic transport component included in the software routes the runoff through a system of pipes, channels, junctions and outfalls at a sub-watershed level. A simulation period consists of several time steps, ranging from seconds to days. For each sub-watershed, the resulting runoff and the flow rate and flow depth for each conduit (pipe or channel) are modelled. Several hydrological and hydraulic processes can be modelled in PCSWMM such as time-varying rainfall, evaporation of surface water, rainfall interception from depression storage, flow routing through conduit networks, surcharging, and surface ponding. SWMM offers 'steady flow', 'kinematic wave' and 'dynamic wave' as routing methods, besides five infiltration methods: Horton, Modified Horton, Green-Ampt, Modified Green-Ampt and Curve Number. To develop the models, we used the 'dynamic wave' routing method, because it solves the complete one-dimensional Saint Venant flow equations and therefore theoretically produces the most accurate results (Rossman and Huber 2015). These equations consist of the continuity and momentum equations for conduits and a volume continuity equation at nodes. Curve Number was used as the infiltration method because it is based on land use classification and hydrologic soil groups. It assumes that the total infiltration capacity of a soil can be found from the soil's tabulated curve number (Cronshey, Roberts, and Miller 1985). During a rain event, this capacity is depleted as a function of cumulative rainfall and remaining capacity (Rossman and Huber 2015).
Additional hydrological modelling was carried out using the software HEC-HMS Version 4.7 (Scharffenberg 2016) for comparison. This simulates the surface runoff response of a watershed to precipitation by representing it with interconnected hydrologic and hydraulic components (Oleyiblo and Zhi Jia 2010). The HEC-HMS basin model comprises three vital processes: rainfall abstractions, rainfall-runoff transformation and base flow, whereby in our case the latter was not considered due to the low Quebrada Seca water flows. Snyder's synthetic unit hydrograph was used to transform excess rainfall into runoff.
The rainfall abstractions for both models were estimated using the National Resources Conservation Service (NRCS) method, formerly known as Soil Conservation Service (SCS) method (Appendix B).
The CN values of each land cover class were obtained from standard runoff curve number tables. For high vegetation, low vegetation, pastures and urban cover, the CN values were estimated to be 55, 73, 68 and 92, respectively, also taking soil types into account. Based on the present soil composition, the hydrologic soil group for the Quebrada Seca watershed is Group B (Oreamuno and Villalobos 2015).
A comparative analysis was performed between the HEC-HMS and PCSWMM models, with a view to validating results based on conditions observed on-site. Both models show the hydrological response for the whole Quebrada Seca watershed and use the same base information on land cover and physical characteristics. However, the PCSWMM model was spatially further discretized to obtain a greater level of detail at different control points along the Quebrada Seca River. Figure A2 in Appendix A shows the hydrograph at the Highway 1 control point for each model as a comparative result.
In the PCSWMM model, the entire area was subdivided into 17 elementary drainage areas (Chen et al. 2021) based on topography and known information from the different elementary stormwater drainage systems ( Figure 2). This allows the identification of specific areas where most runoff is generated over time, enabling solutions to be found that involve the parameters with the highest sensitivity level.
Due to the absence of flow-gauging stations in the river, the analysis was carried out using historical precipitation data from previous research, albeit validated via on-site flood-related information, testimonies and experiences of local people. This is considered a valid approach to decrease the uncertainty of hydrological models in ungauged basins. The results of our models were calibrated based on the runoff hydrographs obtained in the studies of (Chen et al. 2021;Oreamuno and Villalobos 2015). These previous studies were used solely to adjust the model for the last year considered, 2019. For all previous years, model parameters remained unchanged. The only exception was the curve number, which was adjusted for  the land cover present in these previous years, thereby highlighting the effect of land cover change on the watershed's hydrologic response.

Documentation of historic flood events
To demonstrate the Quebrada Seca watershed's vulnerability to flooding, a review of the reported flood events in the period 1945 -2020 was carried out. This was accomplished through a review of national newspapers, in search of news on floods in different sectors of the watershed. Table 3 lists the date, the affected area(s) and a brief description of the flood events reported. It is noteworthy that, prior to 1987, there were no reported flood events. This does not imply that prior to that year there were no floods. What is more likely is that the floods did not affect the population or infrastructure to any great extent, therefore not making them relevant news.

Analysis of land cover change in the watershed
The results of the land cover analysis are shown in Figure 3 (  and Figure 4 (1997-2019). In 1945, according to the land cover classification of the orthophotos, the watershed's pervious surfaces consisted mainly of coffee farms, covering approximately 58% of the total area, while urban land cover accounted for just 1.6% (Figure 3). The remaining area was mainly pasture for raising cattle. In the case of the coffee farms, the foliage of coffee plants is excellent for intercepting precipitation. In addition, within the Quebrada Seca watershed, shade-grown coffee farming (illustrated as high vegetation in Figures 3 and 4) has historically been practiced, further increasing the degree of rainfall interception. Therefore, the watershed's hydrological response in the mid-twentieth century was probably much less than under the current conditions governed by urban development. Urbanization began in the canton of Heredia, home to the city of Heredia, and gradually spread north and west to the peripheral cantons of Barva, Flores, and Belén, connecting with their expanding city center nuclei (Figure 3). The watershed's urbanized areas rose from 4.7% in 1960 (orthophoto) to 10.4% in 1975, with most of the increase in the urban centers of the aforementioned cantons linked by a basic road network.
By 1985, 29.0% of the total area had become impervious. Over the ten-year period between 1975 and 1985, urban areas increased by 4.25 km 2 (18.6% of the watershed). This was the decade of greatest urban growth, especially between 1979 and 1985 when urban areas increased by 3.26 km 2 (17.5% of the watershed) with an average urban growth rate of 2.92% per year over the six-year period. This growth was fostered by improved transport networks between the existing urban centers, favoring urban absorption between Heredia and the rest of the cantons. From 1985 onwards, urbanization, both residential and industrial, continued steadily, with an average annual growth rate of 1.08%, to reach 65.6% in 2019. The last major increase was detected between 1997 and 2007, with 3.63 km 2 covered over in this 10-year period (Figure 4).
The conversion of pervious areas originally used for coffee growing and cattle farming into urban land cover occurred mainly in the middle and lower parts of the watershed. However, during this period of urbanization, the pervious areas themselves underwent changes, with the results of our land cover classification showing that formerly shade-grown coffee cultures were converted to non-shade-grown cultures, while many pastures were transformed into new coffee-growing areas. The latter development was positive from a hydrological perspective, since the foliage of coffee plants is better at intercepting precipitation than pastures. Occurring before urbanization accelerated in 1975 (Figure 3), this type of change took place mainly in the upper parts of the watershed, in the cantons of San Rafael and Barva. These changes in pervious land cover were not however significant at watershed level, as they were cancelled out by the ongoing urbanization. They were more relevant for the local-level hydrological response.

Hydrological response due to land cover changes at watershed and sub-watershed level
Accelerated urban growth in the Quebrada Seca watershed, especially from 1979 onwards, has led to an increase in reported flooding events, mainly in the canton of Belén at the downstream end of the watershed (Table 3). Runoff volumes resulting from the impermeabilization of the watershed tend to be higher in the lower parts of the watershed due to the accumulation of downstream flows.
The Quebrada Seca and Burío River sub-watersheds immediately upstream of the confluence were analyzed independently in order to evaluate their respective hydrological responses to land cover changes. Both feature a similar behavior, with urbanization growing slowly until 1979 (slightly higher in the Burío River subwatershed) and accelerating from 1979 onwards, from 9.7% to 21.6% in the Burío River sub-watershed, and from 2.3% to 11.5% in the Quebrada Seca sub-watershed. This accelerated urbanization between 1979 and 1985 significantly affected runoff volumes and peak flow rates, specifically increasing the runoff volume by about 56,000 m 3 (28%) and the peak flow rate by 5.8 m 3 /s (51%) in the Quebrada Seca sub-watershed. In the Burío River sub-watershed, the greatest increase in runoff volume (29%) and peak flow (60%) occurred earlier, between 1960 and 1975. Accelerated urbanization led to a faster hydrological response of the sub-watersheds, as evidenced by the decrease in time-to-peak, which in 1945 was 75 minutes for the Quebrada Seca sub-watershed and 85 minutes for the Burío River sub-watershed; by 2019, it had dropped to 65 minutes in both sub-watersheds. The changed response can also be seen in the maximum flow generated per square kilometer, i.e. the maximum specific discharge: in 1945 it was 2.0 (m 3 /s)/km 2 for the Quebrada Seca sub-watershed and 3.2 (m 3 /s)/km 2 for the Burío River sub-watershed, whereas in 2019, it had tripled in both cases, reaching 6.7 (m 3 /s)/km 2 for the Quebrada Seca subwatershed and 9.8 (m 3 /s)/km 2 for the Burío River sub-watershed ( Figure 5).
The hydrological response at the control point on Highway 1 featured a behavior similar to the above-discussed subwatersheds. Here again, urbanization developed slowly in the period 1945-1979 before accelerating between 1979 and 1985, with the percentage of urban area increasing from 7.2% to Between 1979 and 1985 urban land cover in the whole watershed increased by 15% points, from 14% in 1979 to 29% in 1985. In this same period, total runoff volume increased by 17% and peak flow rate by 26%. The time-to-peak decreased by 25 minutes from 1945 to 2019, reaching a value of 100 minutes in 2019. The maximum discharge generated per square kilometer again tripled from 2.9 (m 3 /s)/km 2 in 1945 to 9.1 (m 3 /s)/km 2 in 2019 ( Figure 5).

Hydrological response due to the land cover change at the elementary drainage area level
Although our analysis of the hydrological response at subwatershed level due to land cover changes revealed differences between the watershed's upper, middle and lower parts, a more detailed assessment at elementary drainage area level enabled an explanation of the impact of urbanization in an even more spatially explicit manner. It also allowed the recognition of irregularities -in most cases abrupt changes at specific sites and periods which subsequently disappeared, restoring the general trend of the data set. The results per elementary drainage area, i.e. the percentage of urban area, total runoff volume, peak flow rate, and specific discharge for the beginning, middle, and end of the evaluated time period, are summarized in Table 4. A more detailed analysis was performed on six representative elementary drainage areas to describe spatially differing phenomena of urbanization and hydrological response (Figure 6).
The percentage of urban area, runoff volume, peak flow rate, and specific discharge for all elementary drainage areas in 1945, and their respective percentage changes for 1985 and 2019 compared to base year 1945, are presented in Table 4. For 1985, all four parameters analyzed show a continuous increase in most of the areas, with the exception A0, A1 and A2, which show a decrease in runoff volume, peak flow rate and specific discharge, despite their percentage of urban area increasing. Areas A0, A1 and A2 are located in the headwaters of the Burío River sub-watershed. In 1945, their predominant land cover was pasture and low vegetation, with small sectors of high vegetation. By 1985, despite the increase in urban coverage, there had been a significant change in land cover from pasture to high vegetation, increasing interception, waterlogging and infiltration, and thus explaining the decrease in runoff volume, peak flow rate and specific discharge. By 2019, all parameters increased considerably for each of the 17 elementary drainage areas. On average, urban land cover increased by 65% points, runoff volume by 80% and peak flow rate and specific discharge by about 220% compared to the base values for 1945.
A different behavior was found in area A7-2, where an irregular increase in urban land cover from 14.7% to 37.7% occurred between 1979 and 1985, only to decrease to 14.4% by 1991 due to an increase in the use of land designated as pasture land (Figure 4). Another example of irregularities is observed in elementary drainage area A0, where an irregular increase in peak specific discharge was noted in 1975 while the percentage of urban area remained constant; however, there was an increase in the curve number and runoff volume, which subsequently decreased in 1979. This increase was due to the change in land cover from high vegetation to low vegetation in 1975, which by 1979 had returned to a condition similar to that of 1960.
The typical development of most of the elementary drainage areas in the lower-middle part of the watershed is characterized by gradual growth of urban areas from 1945 to 1979, but a sudden increase in 1979. Following this abrupt growth, there are two scenarios: the first shows a continuation of accelerated growth in the following years and the second shows a slow growth that tends to be an almost constant value of urban cover. Elementary drainage area A4-1 features this type of behavior: until 1979, the percentage of urban land cover did not exceed 2%. By 1985, following a period of accelerated urban growth, this percentage had risen to 39%, causing a 53% increase in runoff volume and doubling the maximum specific discharge for that area. Although the chart for area A4-1 ( Figure 6) shows a decrease in the percentage of urban areas from 54. 7% in 1997 to 48.7% in 1997, by 2001 this had risen to 60.5% and continued growing steadily, reaching a value of 65% in 2007 and remaining constant until 2019. The chart for elementary drainage area A11 ( Figure 6) shows, from 1979 onwards, a different pattern, with three periods of accelerated urbanization. The results for the first period from 1979 to 1991 show an increase of urban land cover from 0.42% to 39.4%, for the second period from 2001 to 2007 an increase from 50.7% to 72%, and for the third period from 2013 to 2019 an increase from 72.3% to 84% (Figure 4). Elementary drainage area A6 features increased urban growth from 1979 onwards, with two periods of accelerated growth, the first between 1979 and 1991 when the percentage of urban coverage reached 28.9%, and the second between 1997 and 2007 when the urban area doubled in size to reach 67.8% of the total area. By 2019, this percentage had reached 77%, generating a total runoff volume of 163 million liters for the modelled precipitation event. Elementary drainage area A6-1 features constant growth in urban land cover since 1975, with an increase of close to 10% points observed for each period analyzed, reaching 85% in 2019 ( Figure 6). Together with the second largest elementary drainage area, this makes it the area with the highest runoff volume generation: 194.7 million liters in 2019 (Figure 7). Area A6-1 is located where the cantons Flores, Heredia, and Belén meet at Highway 1, with each of the cantons having a more or less equal share of the drainage area. The area is characterized by the intensive development of industrial areas with large parking lots.
In elementary drainage area A4, urbanization began earlier. The results show that urban coverage was 7. 6% in 1945, rising to 20.5% by 1960, and doubling to 44.8% by 1975. In each of the following three periods, the increase was approximately 10% points in each period, reaching 73.3% in 1991. Growth has since continued, reaching the highest percentage of all elementary drainage areas (88.7% in 2019). In combination with its larger area, this makes it one of the areas with the highest runoff volume generation, 163.9 million liters (Figure 7). Elementary drainage area A7-1 features a similar behavior to area A4 in terms of urban growth, with urban land coverage already reaching 7.6% in 1945 and increasing almost linearly until 1975, when it reached 23.6%. In the following four years, it almost tripled, Table 4. Variation of degree of urban area and modelled hydraulic indicators for two time periods (1945-1985 & 1945-2019)  reaching 60.7% in 1979. This abrupt increase caused its maximum specific discharge to rise from 3.3 m 3 /s/km 2 to 9.7 m 3 /s/ km 2 , increasing its runoff volume by 58%. After 1979, urbanization growth rates were moderate, with even small decreases being observed, for example in 2013, when it dropped from 75.4% to 72.5%. In 2019, a further increase in urban coverage was observed, bringing it up to 83.6%. This type of early urban growth can be attributed to the proximity of the elementary drainage areas to population centers, such as Heredia in the case of A4, and San Antonio de Belén in the case of A7-1. Figure 7 categorizes the elementary drainage areas based on their percentage of urban area (indicated by different textures) and their maximum specific discharge (color coding), highlighting those elementary drainage areas with the highest runoff volume for the years 1945, 1985 and 2019. The area with the highest runoff generation varies: in 1945, area A1 ranked first (115.6 million liters), while in 1985 it was area A4 (160.2 million liters), and in 2019, area A6-1 (194.7 million liters). The areas generating the highest runoff volumes are located in the middlehigh part of the watershed. Figure 7 shows the significant urban development of areas A3 and A5 which, despite not being among the largest areas, are those with the highest specific discharge generation per square kilometer, 20.4 m 3 /s/km 2 and 24.9 m 3 /s/km 2 respectively. Although area A1 was the only area in 2019 with less than 20% urban area, due to its size, it was one of those generating the highest runoff volume.

Discussion of the results in the context of other studies
The results of the land cover classification show a steady but spatially and temporally heterogeneous urbanization of all parts of the watershed except the headwaters. With satellite imagery becoming available, shorter land cover assessment intervals and higher-resolution spatial land cover classification became possible, allowing an even more detailed quantification and understanding of land cover change dynamics within the watershed, especially against the background of rapid urbanization. The study focused on urbanization-induced land cover changes between 1945 and 2019. For future studies, satellite images from the Sentinel-2 missions are available. These provide data with an even higher spatial resolution, thereby enabling a more detailed study.
Changing over time, land cover classification was used as a variable in modelling the watershed's historical hydrological response, with all other parameters remaining constant. This procedure enabled a study focused on the sole impact of land cover changes, ignoring other possible environmental changes, e.g. climatic conditions, river morphology or drainage system development. By using a model structure of elementary drainage areas subdividing the watershed into smaller hydrological units belonging to sub-watersheds, it was possible to develop a spatially distributed modelling and analysis of the hydrological response to spatially heterogeneous land cover changes.
The spatial discretization employed in the historical land cover analysis enabled the identification of those areas that have experienced the greatest urban growth and currently present a higher percentage of imperviousness. These areas correspond to parts of the watershed where urban development probably took place more rapidly than the development of infrastructure able to handle the increase in runoff volumes or flow rates. They thus correspond to critical parts within the watershed that could eventually be prioritized when implementing flood mitigation measures. However, the most urbanized areas are not necessarily those most prone to flooding, but rather the downstream areas. This is because the hydrological response of urbanized areas is mainly characterized by conveyance systems designed to transport rainwater as fast as possible to downstream water bodies. This phenomenon, complemented by an increase in urban coverage, reduces the watershed's time of concentration, with even short intense storms (Ruberto et al. 2006) causing the river to exceed its hydraulic capacity, in turn leading to flooding and landslides in downstream sectors.
In the Quebrada Seca watershed, the population in the lower part of the watershed, in the canton of Belén, has been the most affected by flooding over the last 30 years. In the late 1980s, flooding began to become a visible problem in the watershed, coinciding with the peak of urban growth between 1979 and 1985 mainly seen in the middle part of the watershed, and the subsequent sustained growth. Modernization and urban development grew in line with the so-called 'demographic boom' of the 1970s and 1980s, characterized by the creation of public entities, autonomous and semi-autonomous institutions and state growth and leading to an increase in urban coverage generated by population growth (Rodríguez-Argüello 2010). This phenomenon is also visible in other urban watersheds within the GAM. For example, the María Aguilar River watershed has also experienced accelerated urbanization, tripling the peak flow rate and doubling the runoff volume from 1945 to 2013 (Ramírez-Sandí 2016). This situation resembles that in the Quebrada Seca watershed, where the peak flow rate and runoff volume were 65.7 m 3 /s and 103.5 × 10^4 m 3 respectively in 1945, but 208 m 3 /s and 191.5 × 10^4 m 3 in 2019.
In watersheds characterized by a high level of urbanization, measures to counteract the flood-generating effects of soil impermeabilization are usually limited. This is due to the high costs involved in replacing the existing infrastructure with one that meets current hydrological requirements, or due to a lack of space for developing new infrastructure complementing existing facilities. In such cases, the implementation of green infrastructure as Nature-based Solutions promoting rainwater infiltration and retention (e.g. green roofs, bio-gardens or bioretention areas) are a promising alternative, since they can be adapted to fit in with the existing infrastructure -as a form of retrofitting (Hack and Schröter 2021). In the context of the SEE-URBAN-WATER research project, such solutions have been studied in different scenarios for the Quebrada Seca watershed (Singh, Kumar Sarma, and Hack 2020;Arthur and Hack 2022;Chapa, Pérez, and Hack 2020;Fluhrer, Chapa, and Hack 2021). However, in cases where introducing green infrastructure is not feasible or where its effect is not sufficient to mitigate the associated problems, an early warning system (EWS) could be developed as a management tool. In small and highly urbanized watersheds, such as the Quebrada Seca watershed, any EWS should be based on precipitation thresholds, allowing the population, especially in the lower reaches, to be warned of potential flash floods.
In the particular case of the Quebrada Seca watershed, in addition to the hydrological and hydraulic considerations needing to be taken into account implementing flood mitigation measures, socio-political considerations also play a role. Since the watershed covers six different municipalities, measures must be carried out jointly by the municipalities under a shared vision. This adds further complexity to flood management, as any actions taken by municipalities in the upper parts of the watershed (San Rafael and Barva) and the middle parts of the watershed (Flores and Heredia) will have a direct impact on the canton of Belén. Therefore, it is essential that, prior to implementation, cooperative strategies are first established between the municipalities involved to maximize the overall benefit for the whole watershed and its population.

Discussion of the study's limitations
Any analysis of land cover change has an associated uncertainty due to differences in the spatial resolution of the employed images resulting from different data sources and classification methodologies (manual for orthophotos, supervised for Landsat images). The resolution of Landsat images has progressively improved, from a pixel size of 60 m in 1979 to 30 m in 1985, in our case the period of greatest urban growth. Moreover, for 1945 and 1960 there were no Landsat images yet available, though panchromatic orthophotos were available. While these images had a higher spatial resolution than the Landsat images, their classification had to be performed manually. Nevertheless, we considered it relevant to use these orthophotos as they provide extremely valuable information about the watershed's land coverage prior to mass urbanization. The resolution heterogeneity of the images limits comparisons between certain years, since classifications are a function of spatial resolution. However, this limitation was mitigated by defining a reduced number of land cover classes, thereby maximizing the contrast between them and meaning that classification is less affected by changes in the spatial resolution of the images. Given its hydrological focus, we considered that the four classes defined were sufficient to characterize the runoff response of the Quebrada Seca watershed.
Any historical analysis of land cover changes has the limitation of it not being possible to verify whether classification for past years is appropriate, since no field corroboration can be carried out. The uncertainty associated with this limitation was reduced by comparing our results with other studies focused on historical land cover changes in sectors within or close to the watershed (Pujol-Mesalles and Pérez Molina 2013). Though not having a hydrological focus, such works allowed the historical classification of land cover changes in the Quebrada Seca watershed to be validated. Moreover, at watershed level, the land cover classification in the period analyzed is consistent with previous studies analyzing specific years (Masís-Campos and Vargas-Picado 2014). The representativeness of the land cover classification at watershed level in turn allowed the resulting classification for the sub-watersheds and elementary drainage areas to be validated.
The uncertainty of the model is an important limitation that must be considered in any hydrological modelling, since a model has to use estimates for certain periods that are based on meteorological information constantly evolving due to climate change. A further limitation is the necessity to make assumptions due to a scarcity of information. There was for example a question-mark over the accuracy of the available topographical data used to determine constants and the watershed's morphological parameters.
The lack of detailed information on stormwater drainage systems is a further important limitation that could considerably impact results. While the aim of the study was to analyze the increase in runoff due to land cover changes, there is no guarantee that flooding is solely caused by the increased percentage of impervious land. Factors such as channel modifications, lack of hydraulic capacity at bridge crossings or the nature of the storms caused by specific climatic conditions can also be the cause of flooding and slope instability problems in the watershed.

Generalizations of results
The historical land cover analysis and the hydrological modelling presented in this study largely relied on freely and publicly available data and software, except for the orthophotos and the precipitation data required for the hydrological modelling. The study can thus be replicated in most other places to provide insights on historical land cover changes, not even limited to urbanization processes. In combination with precipitation events of interest, the hydrological response due to such changes can be modelled. Such an analysis can contribute to identifying tipping points in land cover changes and the subsequent hydrological response, thereby guiding, at least with regard to preventing urban flooding, a sustainable urbanization approach. In this respect, the results of this study can directly inform authorities about the impacts of urbanization on the hydrological response of watersheds in other regions with similar conditions and characteristics, for example in several other highly urbanized watersheds of similar size in the same Greater Metropolitan Area of San José in Costa Rica (Pujol-Mesalles and Pérez Molina 2013), but also in other metropolitan areas of Central and South American countries.
The watershed studied in this work covers six municipalities (cantons), each of them formally in charge of urban planning within their administrative boundaries. Shared municipal responsibility is common in urban watersheds, where a spatially detailed analysis of urbanization and its hydrological consequences can support better decision-making to avoid flood-promoting urbanization and achieve better coordinated and sustainable development. A spatial discretization of the hydrological model into drainage areas entirely located within a single municipality can for example enable the municipality concerned to manage the explicit hydrological response through a targeted planning and design of urbanization and stormwater management.

Conclusions
The objectives of this study were a long-term analysis of land cover change as a result of urbanization and an analysis of the flood generation due to land cover changes using a hydrological model. The historical variation in the hydrological response of the Quebrada Seca watershed is mainly dictated by urbanization, which has led to an increase in runoff volumes (+80%), higher peak flow rates (+220%) and shorter times-topeak (−25 min) between 1945 and 2019.
Although changes in land cover in the upper part of the basin between 1945 and 1960, mainly in agricultural use, also resulted in runoff reductions; however, they did not compensate for the general increase in the watershed's hydrological response in magnitude and timing.
Urban growth between 1945 and 1975 was gradual, with the urban area increasing from 1.6% to 10.4%, but from 1975 to 1985, the percentage of urban land use tripled reaching 29%.
The increase in impervious cover was mainly to be seen in the central part of the watershed, where currently the areas generating the highest runoff volumes are located. By 1985, accelerated urban growth had caused a 20% increase in the total runoff volume of the Quebrada Seca watershed compared to 1975, increasing from 1,135,000 m 3 to 1,359,000 m 3 , indicating that a tipping point of urbanization was reached in the late 1980s, resulting in more frequent flooding from that moment on, as also documented in newspaper articles. This study addresses the relationship between the flooding and accelerated urban growth in cities, especially when such growth does not include measures to mitigate increased runoff volumes generated by soil sealing. To evaluate a possible relationship, it is necessary to compare increased surface runoff with the hydraulic capacity of the watercourses in question, with the objective of establishing what percentage of imperviousness begins to represent a dominant factor for flood generation in a watershed's different sectors. This hydraulic approach is considered to be an important contribution to future work in this field.
To avoid flooding problems, urban growth should be accompanied by mitigation measures for hydrological processes that are altered during land cover change. Introducing structures designed to imitate lost (natural) hydrological cycle processes within an urban drainage area, such as green infrastructures, are a possible solution. As such measures normally increase the time of concentration through retention and reduce runoff volumes through infiltration and evapotranspiration, they should be further investigated and considered when making decisions on ways of reducing flooding, but also regarding possible other co-benefits for society and nature.