Hydrodynamic-Based Numerical Assessment of Flood Risk of Tamuín City, Mexico, by Tampaón River: A Forecast Considering Climate Change

: Climate change has unchained several natural extreme phenomena, including a major frequency and intensity of ﬂooding episodes. From these, the ones of greatest importance are those which endanger human settlements as well as socioeconomic activities. This is the case of Tamuín city, settled in the shore of Tampaón River, in Mexico. In this work, we performed a detailed numerical modelling of the hydrodynamics of the zone, considering in situ topographic and bathymetric data as well as hydrodynamic parameters. Severe rainfall scenarios were simulated in order to determine the zones which are prone to ﬂooding, as well as the potential periods of time between the beginning of the rainfall up to the ﬂooding, considering the potential effects of climate change in the precipitation rate. The outcome of this research will help local governments undertake preventive actions to reinforce the identiﬁed risky zones, thus providing an adequate protection of rural and urban zones, as well as their inhabitants and their economical activities from current and future ﬂoods, considering potential climate change effects.


Introduction
Since the 1950s, Global Warming (GW) has proven to be driving a Climate Change (CC) with severe present and future consequences for mankind as well as for nature itself [1]. Global warming is defined as the rise of the global mean temperature of air over land and ocean [1], which affects the whole Earth's ecosystems. Undoubtedly, this temperature rise is the result of the modification of Earth's global energy balance, increasing the total global heat of the planet [2]. CC has unchained a very intense research activity since the second half of 20th century. Furthermore, in the last two decades, a very important branch of such investigation has focused in providing projections or predictions of the global mean temperature behaviour in diverse scenarios. Modelling has had a major role in these predictions. One of the most clearly affected of Earth's systems due to CC is the planet's hydrology. Much research on the impact of CC in large water bodies has been performed, at both the oceanic level [3][4][5][6][7] and for inland surface waters (lakes, lagoons, estuaries, wetlands, etc.) [8][9][10][11][12][13][14][15][16][17][18][19][20][21][22], groundwaters [23][24][25][26][27][28], etc. Of particular interest from the human point of view are rivers, due to from five GCMs, finding increased mean and maximum river discharges in the northwest while decreased ones in eastern basins. Patrick et al. [40] obtained the 100-and 500-year future flood maps as two-dimensional delineation of potential flood extent in New York City, identifying three different zones: areas that shall certainly flood, areas that shall not, and areas that might flood depending on sea level rise projections through the 21st century. Arnell and Gosling [41] assessed the implication of CC for global river flood risk by region, based on the estimation of flood frequency relationships at a grid resolution of 0.5 • × 0.5 • , using a global hydrological model with climate scenarios derived from 21 climate models, together with projections of future population for 2050. The outcome shows a strong regional variability (Asia being more affected), and considerable variability between climate models with a range in increased exposure of 31-450 million people and 59-430 thousand km 2 of cropland, with a change in flood risk varying between −9 and +376 over the risk in 2050 in the absence of CC. Another work performs global future river flood risk projections that separate the impacts of CC and socio-economic development, based on an ensemble of climate model outputs, socio-economic scenarios, and a hydrologic river flood model combined with socio-economic impact models [42], and finds that, globally, absolute damage may increase by up to a factor of 20 by the end of the century without action, with Southeast Asian countries facing a severe increase in flood risk as consequence of CC, while African countries facing a strong increase in risk mainly due to socio-economic change.
The generalised predicted increment in global flooding risk due to CC all around the world has raised the major concern of scientists as well as policy-and decision-makers on how to mitigate the effects of CC in floods of water bodies, particularly of those near human settlements. With respect to the minimisation of the impact of CC in river flooding, Francesch-Huidobro et al. [43] examined the nature of visions, strategies, plans and programmes implemented by the cities of Hong Kong, Guangzhou and Rotterdam, which are highly exposed to flooding because they are located in river deltas. The work by Alfieri et al. [44] is devoted to evaluating adaptation strategies to limit the impact of river flooding on population and assets at European scale, under a high-end global warming scenario over the time range 1976-2100. The results show that expected damage and population affected by river floods can be compensated through different configurations of adaptation measures, which should favour measures targeted at reducing the impacts of floods, rather than trying to avoid them. In [45], the authors evaluated the impact of CC and a high urbanisation rate in Addis Ababa through SWAT, and found 10% and 25% increases in peak flows due to these two factors, respectively. In [46], a tracing of the application at municipal level of the Safety Plan for 100 mm/h-Rainfall, a scheme for preventing and mitigating inundation caused by extremely heavy, short-term rainfall established by the Japanese central government, is performed. Prospects of how a safety plan should be implemented in an urban watershed are discussed. The authors of [47] developed the Coastal City Flood Vulnerability Index (CCFVI), based on exposure, susceptibility and resilience to coastal flooding, to link theoretical flood vulnerability concepts and the day-to-day decision-making process. They then applied it to nine different cities under both current and CC conditions, demonstrating it to be an effective tool to obtain a broad overview of flood vulnerability in these conditions, with regard hydro-geological, socio-economic and politico-administrative components. Jeffers [48] evaluated policy-and decision-making responses to flood hazards heavily based in physical exposure protection and largely neglecting socio-economic vulnerability, and found them to be ineffective in facilitating effective adaptation to floods and loses. Bleck et al. [49] described the adaptation of flood protection with respect to CC and global sea level rise of the State and City of Bremen, in the General Plan for Coastal Protection-Lower Saxony/Bremen. It outlines administrative aspects as well as special technical requirements for harbour operation, recreational use and urban architecture, among others. In [50], the authors reviewed the evidence about the combined impacts of urbanisation and climate on urban flooding and water quality of inland catchments of the United Kingdom, and assessed its utility for setting environmental legislation and managing the urban water environment. Medium-high confidence evidence show an increase in pluvial and fluvial flood risk, and further reduction in water quality caused by point source pollution and altered flow regimes. On the other side, identified factors that limit the utility of such evidence for managing the urban environment are CC projection uncertainty and suitability, lack of sub-daily projections for storm rainfall, complexity in managing and modelling the urban environment, and lack of probable national-scale future urban land-use projections. In addition, Ek et al. [51] presented the different landscapes along Sweden and considered both fluvial and pluvial floods which are extreme but with low probability and have relatively small, but significant in financial terms, risks, respectively. This analysis considers that, as result of CC: floods are expected to occur with increased frequency in Sweden; the expected temperature rise in Scandinavian countries is expected to be greater than the global average; and precipitation is expected to increase mainly in the north and in the southwest, while the southeast part of the country will face increased periods of drought.
The case of study in this work is located in the Huasteca Potosina region, in northeastern Mexico, with a surface of about 11,000 km 2 . This region stands out for the abundance of water and biodiversity resources [52][53][54][55][56][57][58], it possesses mainly tropical rainforest, semi-tropical and sub-humid tropical climates [52,54,59] and it has sparse and small populations dedicated mainly to agricultural and cattle farming [60][61][62]. These rural societies, formed by an important percentage of indigenous people organised in communal lands (ejidos) [63], have undergone in the last two decades very important ecological transformations by human activity due to the implementation of regional development projects [62,64]. The land-cover and land-use changes in the region show a greater fragility of the natural vegetation and the agricultural areas with respect to hydrometeorological phenomena and plagues [64], which threaten local activities. This fragility is expected to be dramatically increased by CC. Nevertheless, mainly due to the lack of large human settlements and of industrial activities, this region has gone unnoticed by scientists and by federal government, thus there is very little knowledge of the current and future impacts of CC in the hydrological processes of the region. A particular case of interest within Huasteca Potosina is Tamuín city, settled at the borders of Tampaón River. This small city, inhabited by about 38,000 dwellers, stands out for being prone to flooding for several centuries, according to historical records. In 2003, an early study by Hudson and Colditz [65] measured that the flood extent in Tamuín valley can reach up to 77.3% of the total surface of the valley (207.9 km). Thus, it is mandatory for a study to be performed to accurately model flood risk considering potential effects of CC in local hydrological processes and be able to undertake the necessary actions to protect the lives and socio-economics of the local population.
In this work, we performed an assessment of the flooding risk of Tamuín city, through a parallel self-made hydrodynamic model previously validated in diverse studies [66][67][68][69][70] and carefully calibrated with field measurements, with the idea of having full control of processes, methods and precision provided by owning the code. Streamflow of Tampaón River in wet and dry seasons was simulated according to historical records. The results coincide with measurements with a Nash-Sutcliffe efficiency of 0.57 [71] and were compared to the predicted streamflow considering the regional effects of CC. To take into account CC, we considered projected precipitations for mid-century (2046-2065) periods obtained with the ensemble of climate models participating in the Coupled Model Intercomparison Project Phase 5 (CMIP5) [72] under the stabilisation global warming (RCP4.5) scenario [73], provided in IPCC AR5 [74]. This article is organised as follows. In Section 2, we present a detailed characterisation of the region under study. In Section 3, we expose the methodology applied, including field measurements and hydrological conditions, as well as the calibration of the model. The results of the flooding of adjacent lands without and with considering of CC effects are presented in 4. In Section 5, the conclusions of this study are presented as well as some final important remarks. Finally, the applied hydrodynamics and the numerical implementation of the model can be consulted in Appendices A and B.

Region Under Study
As previously mentioned, the region under study is the Tamuín municipality, which has suffered serious flooding episodes during about 500 years according to the available records.

Location
The region under study is located in a section of Tampaón River which passes through Tamuín municipality in the state of San Luis Potosí, México. It lies southeast of the Eastern Mexican Belt, and it consists of a coastal plane. The limits of the studied zone are between 21 • 46' and 22 • 24' N latitude, as well as between 98 • 24' and 98 • 27' W longitude (see Figure 1). The zone under study is embedded in the Huasteca region, in San Luis Potosí state, Mexico, with an altitude of 30 masl, and it belongs to RH-26 Pánuco hydrographic basin [75,76]. Important water bodies in this basin are Los Patitos, Tansey, Brasil, San José del Limón, Palmas Cortadas and Mirador lagoons.
Tampaón River crosses Tamuín city, and it later joins the Moctezuma River 4 km downstream after the city, to form Pánuco River, which is one of the five main rivers of Mexico ending at the Gulf of Mexico [65,75].

Evolution of the Region
Tamuín is a human settlement dating from the pre-Columbian time, where the Huastec civilisation flourished due to the welfare that the fertile lands provided by the waters of Tampaón River. In the south (Figure 1), the location of Tamtoc archaeological site can be observed surrounded by a portion of Tampaón River; it dates from 900 A.D., and is estimated to have fallen about 1100 A.D. [77]. Afterwards, the Huastec culture moved a bit to the northeast, and, by the 15th century, they founded Tamohi settlement at the right shore of the river (see Figure 1) [78,79].
At the time of the Spanish conquest, the village of Tamuín was settled about 2 km to the north of Tamohi archaeological site (see Figure 1). Over time, it became clear that this location was not a good choice because the zone is very prone to flooding, destroying economic activities such as agricultural and cattle raising, as well as households and human lives. Although several historic events have changed the political limits of Tamuín, it was not until 1892 that the village was refounded in its actual location (on the north shore of Tampaón River, see Figure 1), in part to diminish the flooding risk.
Nonetheless, in 1955, Tamuín village suffered an intense flooding, which remained for more than three weeks, a consequence of Gladys and Hilda hurricanes, which also destroyed the only road connecting the village. Although the shore was raised, unexpected floods still happen in the town, which have become more frequent and more intense due to CC.
Nowadays, Tamuín has important natural resources and a population of almost 38,000, most of whom are native and Mestee, as well as a multi-faceted economy [80]. Ninety per cent of the local economic activities are tourism, farming and cattle raising. Along water bodies such as Tampaón River, which is under study, fishing and commerce are also important economic activities. These activities have changed the hydrographic features of the rivers, modifying its hydrology and its water quality [81].

Hydrodynamics
Of the two weather stations within RH-26 Pánuco hydrographic basin, data from 24,139 Tamuín DGE station were used because it has been under operation for the last 35 years. The referred station is located at the coordinates −98.81194444444439, 22 To estimate the flow in the reach of the river under consideration (see Figure 3), we used the equivalent TUH (Triangular Unit Hydrograph) method to compute the Curve Number (CN) [82,83]. These methods, developed by the USBR (United States Bureau of Reclamation) and the SCS (Soil Conservation Service), allow estimating flow rates by dripping for different return periods (T r ), which can be observed in Table 1, while a TUH for a return period of 20 years is shown in Figure 2.

Climate
Weather is warm during most part of the year, ranging from cool to cold between November and February. According to Tamuín DGE station, the minimum registered temperature is 17 • C, the maximum is 34.6 • C, and the average is 25.4 • C.
Annual average precipitation according to the station is 1007.5 mm, but for the whole hydrologic basin is 1427.53 mm, according to the Thiessen polygon method [82]. While the dry season begins on 19 April on average, the wet season usually starts on 8 August [60].

Methodology
The methodology herein presented consists in the development of a hydrodynamic numerical model to simulate the hydrodynamical behaviour of the section of Tampaón River under study in the presence of periods of high precipitation, in both current and intensified by CC conditions [84,85]. Along this section, we describe the performed measurement campaigns in order to obtain reliable hydrodynamic data (flow rate, surface elevation and flow speed), as well as the implementation and calibration of the method along with the particular setup of the numerical model for the case of study. The making of the model, from the governing equations up to the numerical solution strategy are discussed in Appendices A and B.

Field Measurement Investigation
In situ measurements were performed on the first day of the wet and dry seasons according to Section 2.4. Both data about the seasonal variation of Tampaón River levels and the hydrodynamic behaviour of the flow, e.g. the flow rate, were acquired.
Hydrodynamic measurements were obtained through an ADCP RiverRay device, which works based on Doppler effect. This device was carefully fastened to a Hydroboard, specifically designed for its stability in water. Data were recorded in real time, along the GPS position where such data was taken.

Transects
To obtain hydrodynamic data, two kinds of transects in the river were performed. First, to measure the flow rate of the river, a control section of the river was selected. Then, the first transect consisted in crossing the river perpendicularly (see Figure 3). On the other side, to obtain the value of the water surface as well as flow velocities, the second transect was designed so to cross the river all across the section under study in zigzag, as shown in Figure 4.

Bathymetry
To obtain the bathymetry, measurement campaign data were processed as follows: surface elevation data obtained in the second transect were used to generate a digital map of geo-referenced isobates. Figure 5 shows these isobates along with level curves of nearby zones with the aid of a topographic survey as well as with previously existent maps.

Flow Velocity
The flow velocity is another hydrodynamic parameter of importance for numerical simulations. To obtain it in the region under study, data from the second transect were processed through WinRiver II software [86]. The obtained flow velocity can be observed in Figure 6.

Los Arenales zone
Turbulence zones

Measurements start
Tamuín -San Vicente bridge m/s As shown in Figure 6, the (measured) riverbed of Tampaón River in such section presents very low speeds (0.05-0.2 m/s), which is consistent with the lower depths of this regions according to the obtained bathymetry. This makes this reach of the river the most secure for navigation while the most insecure for rising river levels, which represents a great threat for the population living above the river's border in this section.
In the same figure, at the other border of the river, higher flow speeds can be observed, with peaks of almost 1 m/s. On the other side, a turbulence zone is identified and located between (22.002083 • N, −98.767435 • ) and (22.001364 • N, −98.765955 • ), which is dangerous for navigation (even for measurement campaigns). Nevertheless, historically, this zone is less prone to overflows. Furthermore, in the lower border of the river, there is no population, because it consists of a zone currently dedicated to agricultural activities (see Figure 6).

Flow Rates
To estimate the flow rate crossing the control section (see Figure 3), data from the first transect were used. For the dry season, a flow rate of 36.39 m 3 /s was obtained, while, for the wet season, it corresponded to 86.08 m 3 /s (see Figure 7).  The set of equations used for the development of the hydrodynamic model, as well as the implemented numerical solution strategies are presented in Appendices A and B, respectively. The reader interested in those details as well as with the software implementation is directed to those appendices.

Boundary Conditions
River's hydrological flow, Q = 86.082 m 3 /s, was considered constant, thus Dirichlet boundary conditions were used (see Figure 8). The outflow boundary condition is of Neumann's type, ∂U/∂η = 0 (see Figure 8). For dry cells, a non-slipping condition was applied, imposing a normal velocity of zero.

Calibration Approach
The calibration phase consisted in the iteration of introducing parameters and initial data to the code for the model to be executed, and comparing the results with those observed/measured in the zone. This process was repeated until the model matched the required error. If there was a considerable discrepancy between modelled and observed data, initial data of the model were modified and the process was iterated. Calibration ended when the difference between modelled and observed data was below the value of the expected error. For this case, the parameters to calibrate were the flow speeds.
For a correct calibration, Chezy coefficients in the field measured speeds were modified until finding their adequate values for this ground, in order to yield to a more precise comparison between modelled and observed speeds. The calibrated Chezy coefficient was 55.30. The friction coefficients were determined by processing satellite images.
To corroborate the quality of the bathymetry, it was compared with a previous one taken by the government agency INEGI [80]. Moreover, the bathymetry obtained for this study was of much higher resolution than the official one. Figure 9 shows the computed flow velocities. These velocities must be compared with the measured ones, shown in Figure 6.
The error between measured and calculated values of flow speed for calibration purposes were performed in the points of the spatial mesh (see Figure 6). For the zone called "Los Arenales", 704 m upstream of the Tamuín-San Vicente bridge, the error corresponds to 14.9% (see Figures 6 and 9). Within the turbulence zone, the lower riverbed, 537 m upstream of Tamuín-San Vicente bridge, the experimental error is 8.09% (see Figure 6). Within the upper riverbed, 246 m upstream of Tamuín-San Vicente bridge, the error is 6.74% (see Figure 6).

Los Arenales zone
Turbulence zones

Measurements start
Tamuín -San Vicente bridge m/s

Future Precipitations in the Zone Under Study
To assess flood risk in the zone for future rainfall events driven by CC, the considered projected precipitations for the next decades are those reported by IPCC [74]. The difference in the 20-year return value of annual precipitation extremes of the 2046-2065 period with respect to the 1986-2005 baseline period for the region under study was 10%. This result was obtained in the ensemble of climate models participating in the Coupled Model Intercomparison Project Phase 5 (CMIP5) [72], for the RCP4.5 scenario, which stabilises radiative forcing at 4.5 Wm −2 in the year 2100 without ever exceeding that value. RCP4.5 includes long-term, global emissions of greenhouse gases, short-lived species, and land-use/land-cover in a global economic framework [73].
To perform comparable simulations between the numerical flood in the zone considering both historical and projected precipitations, the 20-year return values in Table 1 were corrected by considering an increment of 12% in the precipitations, i.e., considering a 2% above the IPCC projected value (10%) as a security factor. Considering this is a small alluvial valley [76], the precipitation values were transformed to peak flows using the TUH method [82].

Numerical Experiments
To better assess the potential influence of CC in the flood risk in the next decades in the region under study, we performed two numerical experiments: the first by feeding the numerical model with historical data of the 20-year return values in Table 1, and the second one considering these values with the mentioned correction, in order to find a critical maximum expected precipitation for the next 50 years.
For the two performed numerical experiments, the initial level of the river was taken as the one measured in the second measurement campaign, as shown in Figure 10, where a 3D view of the surrounding topography can also be observed. The experiments were performed in the wet season to work with extreme flooding events.
The numerical experiments were performed as follows: • A first simulation was set up with historical data with a 20-year return value.

•
The time at which the cross section of Tamuín city settled in the lower part of topography ( Figure 5) was fully flooded was recorded. • A second simulation was set up with the projected streamflows under CC.

•
The time at which the second simulation achieved the flood level of the first experiment was recorded.

•
Comparison between flooding times was performed.
Please note that the flooding times were obtained numerically by integrating the total discharge within the domain.   Historical data Projected data Figure 12. Simulation times for the eight snapshots shown in Figure 11. The blue (circle) line represents the flooding time obtained by simulating the flow rate in an intense rainfall event using historical data (20-year return values), and the red (square) line represents these values intensified by 12% due to potential effects of CC.

Discussion
The identification of the zones prone to flooding, (A)-(G) with respect to time, is fundamental for Tamuín city to perform a critical route map to perform the necessary civil work to safeguard population. This route map is long-term due to the budget limitations of such a small municipality.
As is well known, adaptation plans only based on rising flood protections have the effect of reducing the frequency of small floods and exposing society to less-frequent but catastrophic floods and potentially long recovery processes [44]. In this sense, and with the idea of reducing the frequency of small floods in Zones (A) and (B), the first action undertaken by local government is to erect a gabion-based wall to protect Zone (A), from the beginning of the Choy River mouth up to the beginning of Los Arenales zone (see Figure 13), which is being built at the time of writing this manuscript.

Los Arenales zone
Choy river's mouth Figure 13. First region to be protected by civil work around Zone (A) (see Figure 11).
To create greater resilience to flooding in Zones (B)-(D), public areas are being transformed into flood-resilient urban parks [87]. For extreme flooding episodes, in Zones (F) and (G), which are dedicated to agriculture and cattle farming, the idea is to augment the hydraulic capacity of the main river channels with the installation of meanders, in combination with the installation of dredging stations at the end of these [88].
With respect to the simulated times at which certain flooding stages were obtained for both precipitation scenarios, it can be observed that, in the natural conditions (before any civil work), the south shore of Tamuín city (Figure 11d) gets flooded in approximately 75 min considering an extreme rainfall episode from the historical 20-year return period. Nevertheless, this time is reduced by 23.46% with precipitations taking into account CC projections up to about 57 min.
In the most extreme considered rainfall scenario, the flooding in Figure 11h would occur in about 121 min according to historical precipitation values, while, considering the projections of rainfall intensified by CC, this scenario would be reached in about 86 min, that is, 29.33% faster.
These simulated flooding times are actually being considered by local authorities in order to redesign civil protection plans so they can evacuate inhabitants in enough time in the case of disaster. As for the inhabitants, a plan to induce actions of self-protection is under construction [89].

Conclusions
In this paper, we develop a hydrodynamics-based flooding model taking into account the bathymetry and surrounding topography to assess the flood risk of Tamuín city, Mexico, due to Tampaón River in both historical and climate change projected precipitation scenarios. To achieve realistic hydrodynamics, a horizontal mixing length model for turbulence, as well as shear stresses due to wind and bottom are considered.
Simulations allowed identifying the zones prone to flooding according to different rainfall periods, in both scenarios. Severe flooding episodes are expected to occur between 25% and 30% more quickly considering an increment tendency of 12% for the 20-year return values expected by CMPI5 CC modelling under an RCP4.5 radiative scenario. These results show an important regional increasing trend in the future local precipitations, which contrasts with that recently reported by Ishak and Rahman [90], Ishak et al. [91], who found regional decreasing trends in annual maximum floods in Australia. Results also allow projecting policies to increase Tamuín city resilience to flood, as well as to better draw local civil protection plans.
Finally, it is important to recall that Tamuín municipality is enclosed in a very rural environment, far from large urban settlements. In this sense, this study sets a precedent in the study of the impact of CC in the northeast part of Mexico, in the challenging arena that constitutes the triangle of climatic conditions that join in this zone: the humid mild weather in the east due to the presence of the Gulf of México (in Tamaulipas and Veracruz states), the very dry and desert one in the north (in Nuevo León and Coahuila states), and the steppe weather in the west in the mid zone of San Luis Potosí (in the Eastern Mountain Range, or Sierra Madre Oriental in Spanish). As no previous studies of the effects of CC in this zone exist (to the knowledge of the authors), we expect this work to be a cornerstone to build a robust knowledge of the local impacts of CC so as to be able to protect both nature and human lives, as this geographical region is characterised by a very rich biodiversity, and it is estimated to be inhabited by more than 5 million people.

Future Work
The next step is to widen this study by performing the simulations considering the whole Huasteca Potosina region, which has an approximate area of 32,000 km 2 , to extend the results for them to be useful to a larger population. Another path is to consider long term projections of rainfall for the late century (2080-2100). Acknowledgments: Authors acknowledge EDI grant by Secretaría de Investigación y Posgrado, Instituto Politécnico Nacional, as well as Facultad de Ingeniería, Universidad Autónoma de San Luis Potosí for the facilities to develop this investigation.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A. Hydrodynamics
Appendix A.1. Hydrodynamic Equations The system equations used to calculate the velocity fields are the Reynolds-averaged Navier-Stokes equations for shallow-water flows [70], which are: where U and V are the depth-averaged velocity components in x and y directions, respectively (m/s); g is the gravity's acceleration (m/s 2 ); η represents the free surface elevation (m); ν TH is the horizontal eddy viscosity (m 2 /s); H is the water depth (m); τ w x and τ w y are the wind shear stress terms in x and y directions, respectively; and τ b x and τ b y are the corresponding x and y direction bottom shear stresses. Shear stresses units are m 2 /s 2 .
To compute the free surface elevation (η), we integrate the continuity equation over the water depth. After applying a kinematic condition at the free surface, we have To obtain realistic modelling, a model of turbulence has been introduced to simulate mechanical dispersion phenomenon. A mixing-length model for turbulence is applied [70], which contributes to the horizontal eddy viscosity as: where l h = βl v is the horizontal mixing length (m), β is a dimensionless constant, and l v (m) is defined as: where z b is bathymetry (m), and α = λη κ , λ and κ are dimensionless constants (the later is the von Kármán constant). It must be remarked that this turbulence model was the one that best fitted the measured values of the four mixing-length models that were tested for this study.
On the other side, wind shear stress terms in Equations (A1) and (A2) are given by: where ρ w is air density (kg/m 3 ); ω x and ω y are the wind speeds measured 10 m above the ground level in x and y directions, respectively; and C w is the wind drag coefficient given by the Garrat formulation C w = (0.75 + 0.067ω 10 )/1000, where ω 10 is the magnitude of the wind speed 10 m above the ground level (m/s). Finally, the bottom shear stress terms in x and y directions are computed as: where Cz is the Chezy coefficient. This hydrodynamic model has been calibrated according to Casulli and Cheng [92], and has been tested in [67,68,70,93].

Appendix B. Numerical Solution Strategy
To solve the coupled PDEs system posed in the previous section, a second-order finite difference formulation in both time and space is used. The solution scheme adapts the semi-implicit Eulerian-Lagrangian scheme by Casulli and Cheng [92], treating the advection and diffusion terms differently. While the solution of the non-linear advection terms of Equations (A1) and (A2) is found with a Lagrangian formulation through the characteristics method, diffusion terms are treated in a Eulerian manner through the Adams-Bashforth scheme.
A two-dimensional staggered mesh for the numerical computations is used (Figures A1 and A2). The centre of the cell (black dot) represents free surface elevation (ϕ i,j ), while circles on the faces of the cell represent the velocity components U and V. As the model's spatial resolution has a direct impact in assessing the impact of CC on flooding [94], we have set the mesh resolution to 5 m.  In Figure A1, H i+ 1 2 ,j is the water depth at (i + 1 2 , j), where the U velocity component is evaluated, and H i,j+ 1 2 is the water depth at point (i, j + 1 2 ), where the V velocity component is evaluated. The solution of the system of Equations (A1), (A2) and (A3) is given through a linear systems as follows: U n+1 i+1/2,j = FU n i+1/2,j − g ∆t ∆x η n+1 i+1,j − η n+1 i,j where FU and FV join the advective and the turbulent diffusion terms as: where U n and V n are the explicit velocity components, computed at time n using the characteristics method as: U n i+1/2,j = U n i+1/2−a,j−b = (1 − p) (1 − q) U n i+1/2−l,j−m + qU n i+1/2−l,j−m−1 +p (1 − q) U n i+1/2−l−1,j−m + qU n i+1/2−l−1,j−m−1 (A15) where a and b are the Courant numbers in x and y directions, respectively, from which l and m are the integer parts, and p and q the decimal parts. The present numerical model has been validated and applied in different studies, each of which has consisted in controlled experiments and real free surface flows in rivers, lakes, lagoons, estuaries and nearshore waters. The numerical model was initially validated by León et al. [66] to study the dispersion of saline plumes in coastal zones. It was also used to predict the thermal plume dispersion of the Laguna Verde Nuclear Power Plant, located in Mexico [67]. Barrios-Piña et al. [68] used the code to probe a new turbulence model for flow through vegetation. On the other side, Torres-Bejarano et al. [69] recently used it to assess pollution distribution in an estuarine ecosystem. Turbulence models attached to the hydrodynamic model were tested by Rodriguez-Cuevas et al. [70], where several turbulence models were coded to predict wakes around conic islands, and the results showed good agreement with the measurements experiments.

Software
For the sake of efficiency, the algorithm was coded in Fortran R , and it was parallelised through the OpenMP paradigm. The model (bathymetry and topography) were an input of the Fortran code and they were coded as CSV text files. The succeeding results were visualised with TecPlot R .