Reliability analysis of beaches as defenses against storm impacts under present and future climate change scenarios

Treball realitzat per: Blanca Bevià Brey Dirigit per: Francesc Xavier Gironella i Cobos Vicente Gracia Garcia Grau en: Enginyeria Civil Barcelona, 28 de Gener del 2025. Departament d’Enginyeria Civil i Ambiental TREBALL FINAL DE GRAU Reliability Analysis of Beaches as Defences Against Storm Impacts Under Present and Future Climate Change Scenarios 2 Acknowledgements I would like to express my sincere gratitude to Professors Francesc Xavier Gironella i Cobos and Vicente Gracia Garcia for their guidance throughout this project. They provided exceptional support from the beginning, offering clear direction, necessary resources, and valuable mentorship that were crucial to completing this work successfully. I am deeply thankful to my family for their constant support. Their interest in my research and their commitment in my academic journey have been essential. Their encouragement has helped me overcome challenges throughout this process. My experience studying Civil Engineering at UPC has been transformative. During my time here, I have discovered the fascinating aspects of civil engineering, including its various specializations and practical applications. The bachelor has provided me with both technical knowledge and important perceptions about the field's impact. I have had the opportunity to study with excellent classmates who have become valued colleagues. Together, we have spent many hours working on projects, engaging in technical discussions, and supporting each other. These mutual aids have strengthened my understanding of teamwork and cross-disciplinary cooperation, creating a strong foundation for my future career in engineering. 3 Index Acknowledgements .................................................................................................. 2 Abstract ................................................................................................................... 5 1. Introduction ................................................................................................................. 7 1.1. Motivation ............................................................................................................. 9 1.2. Objectives ............................................................................................................ 10 1.3. Structure ............................................................................................................. 12 2. Theoretical Framework .............................................................................................. 12 2.1. The Ebro Delta: Coastal Characteristics and Environmental Challenges ............ 13 2.2. Coastal flooding on beaches: review of existing equations of wave run-up and flooding on beaches ....................................................................................................... 17 2.2.1. Hunt (1959) ............................................................................................... 18 2.2.2. Yamamoto (1972) ..................................................................................... 18 2.2.3. Holman (1986) .......................................................................................... 19 2.2.4. Kobayashi and Tamura (1986) .................................................................. 19 2.2.5. Mase (1989) .............................................................................................. 20 2.2.6. Nielsen & Hanslow (1991) ........................................................................ 20 2.2.7. Van der Meer & Stam (1992) .................................................................... 20 2.2.8. Stockdon et al. (2006) ............................................................................... 21 2.2.9. Flick and Murphey (2008) ......................................................................... 21 2.2.10. Mase and Shibayama (2009) .................................................................. 22 2.2.11. Summary of wave run-up equations ...................................................... 22 2.3. Selecting the Appropriate Wave Run-up Equation for the Ebro Delta ............... 23 2.3.1.The Iribarren number ................................................................................ 24 2.4. Sea Level Rise Projections and Regional Impact Analysis ................................... 26 2.4.1 Global Context and Mediterranean Specificity ......................................... 26 2.4.2. Regional Downscaling and Mediterranean Projections ........................... 26 2.4.3. Temporal Analysis of Projected Sea Level Rise ......................................... 26 2.4.3.1. Sea levels projections ...................................................................... 27 2.4.4. Comparative Analysis of Projection Sources ............................................ 27 2.4.5. Sea Level Rise Implications and Uncertainties in Coastal Management .. 27 4 2.5. Failure tree: Hierarchical Analysis of Beach Flooding Components and Risk Assessment ..................................................................................................................... 29 3. Code Methodology .................................................................................................... 32 3.1. Introduction ........................................................................................................ 32 3.1.1. Monte Carlo Simulation Model for Coastal Flooding Risk Assessment ..... 32 3.1.2. Astronomical Tide Modelling ..................................................................... 33 3.1.3. Meteorological residue............................................................................... 35 3.1.4. Wave Height Analysis and Weibull Distribution ......................................... 37 3.1.5. Directional Wave Analysis under Storm Conditions .................................. 40 3.1.6. Wave Run-up Calculation (Stockdon et al., 2006) ...................................... 42 3.1.7. Sea Level Rise Scenarios ............................................................................. 45 3.1.8. Total Water Level and Failure Analysis ....................................................... 47 3.1.9. Convergence Analysis of Monte Carlo Simulation Results ......................... 48 3.1.10. Reliability Analysis for Storm-Induced Coastal Flooding .......................... 48 3.1.11. Application of the model: The Ebro Delta System ................................... 52 3.1.12. Beach Characteristics: Site-Specific Analysis for The Marquesa Beach .. 54 3.1.13. Beach Characteristics: Site-Specific Analysis for The Trabucador Beach ...................................................................................................................58 3.1.14. Methodology for Monte Carlo Simulation Variable Conservation .......... 61 3.2.Theoretical Constraints and Assumptions ........................................................... 63 4. Results ........................................................................................................................ 64 4.1. The Marquesa ..................................................................................................... 65 4.1.1. Wave analysis ............................................................................................. 65 4.1.2. Water Level Components ........................................................................... 66 4.1.3. Failure Probability during Storm Events ..................................................... 67 4.1.4. Water Level Exceedance ............................................................................. 68 4.1.5. Convergence Analysis ................................................................................. 68 4.1.6. Reliability Analysis ...................................................................................... 74 4.2. The Trabucador ................................................................................................... 77 4.2.1. Wave Analysis ............................................................................................. 77 4.2.2. Water Level Components ........................................................................... 78 4.2.3. Failure Probability During Storm Events .................................................... 79 5 4.2.4. Water Level Exceedance ............................................................................. 80 4.2.5. Convergence Analysis ................................................................................. 81 4.2.6 Reliability Analysis ....................................................................................... 86 5. Conclusions ................................................................................................................ 88 6. Sustainability analysis and ethical implications ....................................................... 90 7. Sources ....................................................................................................................... 93 Appendices ..................................................................................................................... 97 1. Model code ............................................................................................................ 97 2. Matlab used functions ......................................................................................... 103 6 Abstract This research evaluates the reliability of natural coastal protection systems against flooding under climate change conditions. The study focuses on two morphologically distinct beaches in Spain's Ebro Delta: The Marquesa and The Trabucador. It assesses their protective capacity and reliability during storm events by analysing how their unique geometric characteristics influence their effectiveness as flood barriers. The methodology implements a probabilistic framework using Monte Carlo simulation. This approach integrates astronomical tides, meteorological conditions, and wave dynamics to assess storm impact scenarios. The analysis spans multiple time horizons (2025, 2050, and 2100) across various climate change scenarios (Current, RCP4.5, and RCP8.5), and includes the effect of land subsidence in the Ebro Delta, evaluating both immediate and long-term coastal protection capabilities. The comparative analysis reveals distinct sea-level rise response patterns between the sites. The Marquesa beach currently contributes to effective coastal protection at higher elevations (2.2m), though this reliability diminishes under projected climate change scenarios. The Trabucador demonstrates greater sensitivity to environmental conditions, with its protective function (2.1m) showing accelerated deterioration. Both locations exhibit significant vulnerability in their lower sections (1.3m and 1.2m respectively), with analysis indicating substantial degradation of protective capacity by 2100, particularly under RCP8.5 scenarios. The research provides crucial quantitative evidence for coastal defence planning by emphasizing that preserving beach geometric characteristics through maintenance programs is essential for sustained protective capacity, action particularly critical as climate change continues to reshape coastlines and intensify storm impacts. Keywords: coastal flooding, beach reliability, climate change, sea level rise, storm impact, Monte Carlo simulation, Ebro Delta, coastal protection. 7 1. Introduction Climate change is profoundly transforming our world, irreversibly altering ecosystems, flora, and fauna. One of its most devastating effects is the rise in sea levels, which is already causing more frequent and severe flooding in coastal areas. Combined with the intensification of storms and hurricanes, these phenomena are threatening both communities and coastal biodiversity, leading to an unprecedented environmental and social crisis. As global temperatures rise and weather patterns shift, coastal communities are facing a significant battle against the sea. Coastal flooding events are becoming more common, as storms grow in frequency and intensity. This imminent crisis threatens lives, jobs, and the existence of many coastal communities, the danger of inundation is imminent, and the time for action is now. It is necessary to develop strategies to strengthen the natural defences that have long been protected in the past from the forces of the sea. Beaches play a crucial role in this defence by absorbing the energy of crashing waves and protecting inland areas from the full impact of storm surges. Their effectiveness depends on two key factors: beach width and elevation. When these are compromised, waves can easily surpass the shoreline, leading to destructive coastal flooding. However, as climate change progressively erodes these natural barriers, the question arises: how much protection can the beaches truly provide in the face of future storms? The Ebro Delta, located on Spain's Mediterranean coast, stands as a region of immense ecological, economic, and social importance where the Ebro River meets the sea (Figure 1). It is home to extensive agricultural lands, wetlands, and lagoons that are vital to the local economy and ecosystem. Figure 1: Map of the Ebro River Basin in Spain. Source: Locken, The challenge of managing the Ebro River Basin: Larger than half the countries in the EU. iAgua, https://www.iagua.es. 8 However, the delta exemplifies the vulnerability of coastal areas, as flooding could inundate these regions, causing severe damage to public infrastructure and habitats. The beaches here dissipate wave energy and protecting inland areas from flooding. Yet, the exceptional challenges faced by the Ebro Delta make the role of these beaches even more critical. With rising sea levels and increasingly intense storms, the region's reliance on these natural defences is rapidly reducing, leaving its communities and ecosystems more vulnerable than ever. This study aims to conduct a reliability analysis of sandy beach archetypes along the Ebro Delta coastline, evaluating their effectiveness as natural defences against storm impacts. The analysis will use Monte Carlo simulations to account for the variability in key hydrodynamic and morphological parameters. By modelling a wide range of potential scenarios, the study seeks to assess the likelihood of beach failure and the associated risk of coastal flooding under both current and future climate conditions. This stochastic approach, powered by Monte Carlo simulations, offers a more realistic representation of the variability and uncertainty inherent in coastal processes, in contrast to traditional deterministic methods. The significance of this research becomes even more critical in the context of climate change. As sea levels rise, storm patterns shift, and storm surges grow more intense, the pressure on coastal defences is expected to increase. This, in turn, accelerates beach erosion and diminishes their ability to protect inland areas. As these natural defences weaken, the risk of severe flooding and its economic and humanitarian impacts will intensify. Through this study, the aim is to present fundamental insights into the longterm reliability of beaches as natural barriers. The results of this research will provide essential data to support efforts in coastal protection and disaster risk reduction. As climate change continues to reshape coastlines, ensuring the resilience of natural defences, such as beaches, is becoming ever more crucial to protect coastal communities, infrastructure, and economies. Therefore, this study represents a significant advancement in understanding and addressing these challenges, offering interpretations that can help shape effective strategies to increase coastal resilience in the face of a changing climate. 9 1.1. Motivation Coastal flooding poses an unprecedented threat that demands our immediate attention and action. Around the world, communities, businesses, and vital ecosystems face rising risks as strengthening seas and extreme weather events intensify. The urgency of this challenge, accelerated by climate change, calls for innovative solutions based on rigorous scientific understanding. The research attempts this critical need by investigating how beaches can protect our coastlines. By expanding our knowledge of these natural barriers and their response to climate change, we aim to unlock new possibilities for coastal protection. This work is not just about understanding beaches; it is about safeguarding communities, preserving economies, and protecting the delicate balance of coastal ecosystems for future generations. The Ebro Delta serves as both a scientifically optimal case study and an area of personal significance. Growing up in the region, my family's frequent excursions to the Delta fostered a deep appreciation for this unique coastal environment. These childhood experiences of observing its dynamic landscapes, diverse ecosystems, and cultural importance have developed into a professional commitment to understanding and protecting this vital natural system. This personal connection provides valuable comprehension of the Delta's significance beyond its ecological and economic value. The Delta's characteristics offer a requiring context for investigating coastal protection strategies, with findings that hold potential applications for coastal communities worldwide facing similar challenges. The urgent need to protect our coastlines calls for innovative solutions, and this research rises to that challenge. At a time when coastal communities face climate threats, this research provides vital insights for policymakers and engineers alike. The findings catalyse action, empowering decision-makers to develop robust, nature-based solutions that strengthen our coastlines against climate change. Through this work, we are not just studying coastal resilience, we are actively shaping a more sustainable and secure future for our coastal regions. 16 The Ebro Delta's experience offers valuable insights for developing adaptive management strategies in similar coastal environments worldwide, illustrating the complex decisions required to protect and preserve vulnerable coastal ecosystems in an era of rapid environmental change. The present research focuses on two distinctive coastal areas within the Ebro Delta that exemplify the varying dynamics and vulnerabilities of this complex ecosystem. Marquesa Beach, situated in the northern hemidelta, and Trabucador Beach (Figure 6) in the southern section serve as representative cases for analysing coastal response to climate change impacts. Figure 6: The Trabucador Beach: A Fragile Point in the Time-Critical Battle at the Ebro Delta. Source: RTVE. (2023, May 24). Platja dthe Trabucador: la lluita a contratemps per salvar el Delta de l'Ebre. RTVE. https://www.rtve.es/television/20230524/platja-trabucador-lluita-a-contratemps-deltaebre/2446163.shtml 17 2.2. Coastal flooding on beaches: review of existing equations of wave run-up and flooding on beaches Wave run-up quantifies the maximum vertical elevation a wave achieves above the still water level when interacting with a coastal slope or structural interface. This precise measurement provides comprehensive evaluation of marine energy interaction with terrestrial environments. Hence, by capturing the vertical displacement of wave motion, as presented in Figure 7, run-up analysis offers essential data for understanding coastal vulnerability. The measurement serves as a key diagnostic tool for evaluating the dynamic forces that shape coastal landscape and flood risk. Figure 7: Definition of Wave Run-up. Source: Deborah Villarroel-Lamb (2022). In Quantifying wave runup in data-sparse locations for planning. Journal of Coastal Research. https://www.researchgate.net/publication/359817802_Quantifying_Wave_Runup_in_DataSparse_Locations_for_Planning Reviewing existing equations for wave run-up requires analysing diverse mathematical formulations used to estimate the height that waves will reach on a beach or a coastal structure. These equations are crucial for predicting coastal flooding and designing effective coastal defences. The available formulations vary in complexity and applicability, ranging from empirical formulas for idealized beaches and structures to advanced numerical models for more complex coastal environments. The choice of equation or model depends on the specific conditions of the site and the level of detail required for the analysis. For accurate predictions, it is crucial to consider beach slope, wave characteristics, and potential tidal influences. They aim to provide a comprehensive understanding of wave behaviour and its impact on coastal areas. 18 The lower figure 8 presents a summary of the most used run-up equations: Figure 8: Evolution of the Wave Run-up Equations. Source: Main. 2.2.1. Hunt (1959) One of the oldest run-up equations, Hunt's formula is based on linear wave theory and is typically used for steeper slopes, such as seawalls or revetments. It focuses on wave run-up on impermeable structures or steep coasts. It is simple, and particularly effective for steeper, impermeable coastal structures. Best suited for: Seawalls, steep beaches, rocky shores. [1] 𝑅=𝛾𝐻𝑠tan(𝜃) (Eq. 1) Description: - R = run-up height. - γ = empirical constant (usually around 1.5 for smooth, impermeable surfaces). - 𝜃 = slope angle. - Hs = significant wave height. 2.2.2. Yamamoto (1972) This equation is often used in coastal engineering to estimate wave run-up in various beach conditions, including those with significant wave heights and varying wavelengths, often applicable for preliminary assessments. [2] 𝑅=𝜉⋅√𝐻𝑠⋅𝐿 (Eq. 2) 19 Description: - R = run-up height. - Hs = significant wave height. - L = deep water wavelength. - 𝜉 = Iribarren number (dimensionless parameter). 2.2.3. Holman (1986) This is a simpler equation that predicts the maximum run-up height based on the wave steepness and beach slope. Best suited for: Coastal sites with mild slopes and sandy beaches. [3] 𝑅2% =𝛼𝛽𝑓√𝐻𝑠𝐿0 (Eq. 3) Description: - R2% is the maximum run-up height. - 𝛼 is an empirical constant that depends on the beach slope. - βf = foreshore slope. - Hs = significant wave height. - L0 = deep water wavelength. 2.2.4. Kobayashi and Tamura (1986) Relates significant wave height to wavelength for estimating wave run-up. Primarily used for sandy beaches, this equation is useful when dealing with run-up under conditions of varying wave height and length. It is particularly applicable in coastal areas with consistent wave patterns. [4] 𝑅2% =𝐶⋅𝐻𝑠⋅(𝐿 𝐻𝑠)𝑛 (Eq. 4) Description: - R2% = run-up height. - Hs = significant wave height. - L = deep water wavelength. - C = an empirical constant, which depends on beach slope and wave conditions. - n= is another empirical exponent that accounts for the effects of the beach slope and wave breaking characteristics. 20 2.2.5. Mase (1989) Mase developed an empirical formula based on field data from natural beaches. It is often used in coastal hazard assessments. Uses field data to validate predictions. Best suited for: Natural beaches, mild to moderate slopes. [5] 𝑅 =0,56𝐻𝑠(𝐿0 𝐻𝑠)2𝛽𝑓 0,5 (Eq. 5) Description: - R = run-up height. - Hs = significant wave height. - L0 = deep water wavelength. - βf = beach slope. 2.2.6. Nielsen & Hanslow (1991) This equation was designed for sandy beaches and focuses on storm wave run-up. It is commonly used for beach erosion studies during high-energy wave conditions. Simple and useful for storm conditions. Best suited for: Sandy, erodible coastlines during storms. [6] 𝑅 2%=2,3⋅𝐻𝑠 (Eq. 6) Description: - R = run-up height. - Hs = significant wave height. 2.2.7. Van der Meer & Stam (1992) This formula is specifically designed for rocky slopes and rubble-mound breakwaters. It incorporates the roughness and permeability of the structure, making it suitable for engineered coastal defences. Tailored to rough and permeable structures, such as breakwaters. Best suited for: Coastal defences like breakwaters and revetments. [7] 𝑅 2%=𝜃𝛾𝑠𝐻𝑠𝜀𝑚 (Eq. 7) 21 Description: - R2% = 2% exceedance run-up height. - 𝜃 = slope angle. - γs = reduction factor (accounts for roughness and permeability). - Hs = significant wave height. - ξm = surf similarity parameter. 2.2.8. Stockdon et al. (2006) This model includes both swash and infragravity wave contributions. Useful for predicting extreme run-up. Best suited for: Sandy beaches, low-gradient shores. [8] 𝑅2% =1,1(0,35𝛽𝑓√𝐻𝑠𝐿0+(𝐻𝑠𝐿0(0,563𝛽𝑓2+0,004))1/2 2) (Eq. 8) Description: - R2% = the 2% exceedance run-up height. - βf = the foreshore beach slope. - Hs = the significant wave height at breaking. - L0=gTp 2 2π is the deep water wavelength (g is gravitational acceleration, and Tp is the wave peak period), calculated from Hs. 2.2.9. Flick and Murphey (2008) This model is particularly useful for varying beach profiles where slope and sediment size have significant influences on wave run-up. It is relevant in environments where complex interactions between wave dynamics and beach morphology occur, useful such as in design and safety assessments. [9] 𝑅 𝑚𝑎𝑥 =1,1⋅ 𝐻𝑠 √2⋅𝑐𝑜𝑠𝜃 (Eq. 9) Description: - Rmax = run-up height. - Hs = significant wave height. - θ = slope angle. 22 2.2.10. Mase and Shibayama (2009) This equation is particularly useful for sandy beaches and is focused on empirical relationships, making it relevant in environments where the wave height and wavelength are primary factors. So applied in coastal management scenarios where empirical data is available for calibration, is also useful in designing coastal defences. [10] 𝑅2% =1,5⋅𝐻𝑠⋅(𝐿 𝐻𝑠)0,5 (Eq. 10) Description: - R2% = the 2% exceedance run-up height. - Hs = significant wave height. 2.2.11. Summary of wave run-up equations The development of wave run-up equations over time reflects a progression towards more comprehensive and accurate predictions of coastal flooding. This chronological evolution demonstrates how incorporating a broader range of variables, refining empirical relationships, and advances in the theoretical understanding of wave-coastal interactions have significantly enhanced the precision of these models. Early equations, such as Hunt (1959), focused primarily on key factors like wave height and beach slope. As the field advanced, subsequent equations began to incorporate additional parameters such as wavelength, wave period, and sediment characteristics. Newer equations, including Mase (1989), Flick and Murphey (2008), and Mase and Shibayama (2009), further refined predictions by accounting for varying beach profiles, slope angles, and empirical relationships specific to different coastal environments. This progression underscores the increasing complexity of these models, with each improvement representing a deeper understanding of the dynamic interactions between waves and coastal systems. Including more variables and improving empirical models have led to equations better suited for specific coastal settings, such as sandy beaches, rocky slopes, or low-gradient shores. In conclusion, the evolution of wave run-up equations has resulted in a set of tools that enable more effective coastal management and mitigation strategies. However, the choice of equation depends on several factors, including the specific coastal environment, data availability, and the level of complexity required for the analysis. Some equations, like Holman (1986), are simpler and require fewer inputs, while others, such as Stockdon et al. (2006), are more complex and need a wider range of data. 23 2.3. Selecting the Appropriate Wave Run-up Equation for the Ebro Delta When studying wave run-up for a specific coastal area like the Ebro Delta, it is crucial to choose an equation that aligns with the region's morphological and physical characteristics. Ebro Delta features a low-lying, gently sloping profile with sandy beaches and wetlands, indicating that equations tailored for low-gradient, sedimentary environments are most suitable. While simpler equations like Holman (1986) can provide quick estimates of maximum run-up height based on wave steepness and beach slope, the Stockdon et al. (2006) equation is generally more appropriate for environments like the Ebro Delta. The Stockdon et al. (2006) [8] equation offers several advantages that make it the preferred choice for the coastline being studied: 1. Comprehensive Approach, as the equation includes both swash and infragravity terms, providing a more accurate representation of the processes occurring on low-gradient, sedimentary coastlines. 2. Parameterization, as the equation is parameterized with coefficients that adapt to different beach profiles and wave conditions, allowing for broader applicability across various coastal environments. 3. Adaptability to Beach Conditions, as the Stockdon equation adapts to different coastal conditions by using either the full or simplified version based on the Iribarren number (ξ). The Iribarren number is a dimensionless parameter that represents the type of beach, with low values (ξ < 3) indicating dissipative conditions and high values (ξ ≥ 3) indicating reflective conditions, for dissipative beaches like the Ebro Delta, the full equation is used, while for reflective beaches, the simplified equation neglecting the infragravity term is applied. 4. Suitability for Dissipative Conditions, as the Ebro Delta is characterized by extremely dissipative conditions due to its gentle slope, wide surf zone, fine sediments, and dominance of low-frequency waves. The Stockdon equation, particularly the full version, is effective in capturing run-up dynamics in such dissipative environments. 5. Accuracy: Although all empirical equations have associated uncertainties, the Stockdon equation has a relatively low error of approximately ±0.32 metres for the 2% exceedance run-up height prediction when using the full equation for dissipative conditions. [8] This makes it a reliable choice for the Ebro Delta given the low-gradient, sedimentary nature of the Ebro Delta’s coastline and the equation's strengths in representing swash and infragravity processes, the Stockdon et al. (2006) equation is the most appropriate choice for studying wave run-up in this region. 24 2.3.1. The Iribarren number However, this expression has a variation for extremely dissipative conditions, due to for reflective beaches, the equation for run-up at 2% exceedance elevations can be simplified by neglecting the infragravity contribution (the 0.004 term). Hence, The Stockdon et al. (2006) equation adapts to different coastal conditions by using either the full (Eq. 8) or simplified (Eq. 12) version based on the Iribarren number (ξ) (Eq. 11). The Iribarren number is a dimensionless parameter that represents the type of beach, with low values indicating dissipative conditions and high values indicating reflective conditions. The Iribarren number is calculated as follows (Eq. 11): ξ = 𝛽𝑓 √𝐻𝑠 𝐿0 (Eq. 11) Where: - Βf is the beach slope. - Hs is the significant wave height at breaking. - L0=gTp 2 2π is the deep water wavelength (g is gravitational acceleration, and T is the wave peak period), calculated from Hs. Equations and associated errors based on the Iribarren number: 1. For dissipative conditions (ξ < 3): Full equation: 𝑅2% =1,1(0,35𝛽𝑓√𝐻𝑠𝐿0+(𝐻𝑠𝐿0(0,563𝛽𝑓2+0,004))1/2 2) (Eq. 8) - Error: ± 0.32 metres [8]. Where: - R2% is the 2% exceedance run-up height (the level exceeded by 2% of the waves). - βf is the foreshore beach slope. - Hs is the significant wave height at breaking. - L0 is the deep water wavelength 25 2. For reflective conditions (ξ ≥ 3): Simplified equation: 𝑅2% =0,73𝛽𝑓√𝐻𝑠𝐿0 (Eq. 12) - Error: ± 0.47 metres [8]. It is important to note that the error could vary slightly depending on local conditions, but this provides a general idea of the accuracy. If used in environments like the Ebro Delta, where the coastal profile might not differ drastically from the data set Stockdon equation (gently sloping sandy beaches), the error should be within a similar range. 32 3. Code Methodology 3.1. Introduction A reliability analysis for beach defence systems against storm impacts has been developed in MATLAB, Appendix 1 contains the complete code used, while Appendix 2 provides an overview of the functions utilized within the code. This engineering tool evaluates beach performance under present and future climate scenarios through customizable parameters. The code implements local characteristics such as beach slope, critical water level height, valid wave directions, and subsidence rates to assess sites' reliability and develop evidence-based adaptation strategies. The probabilistic methodology ensures robust reliability estimates through a Monte Carlo approach with a high number of repetitions, accounting for uncertainties in both climate projections (RCP4.5 and RCP8.5 scenarios), subsidence rates, and environmental conditions (astronomical tide, meteorological residue, and wave parameters) from 2025 to 2100. The framework's modular structure enables the evaluation of individual beaches within delta regions by adjusting site-specific parameters. 3.1.1. Monte Carlo Simulation Model for Coastal Flooding Risk Assessment The analysis employs a Monte Carlo simulation model with 1000000 (1M) iterations reaching robust convergence to assess coastal flooding risks. Through random sampling and statistical modelling, the model generates comprehensive estimations of flooding scenarios, providing essential insights for coastal management. The method incorporates uncertainty by assigning probability distributions to input parameters, based on historical data, and statistical analysis. The model evaluates flooding risks and infrastructure impacts while accurately representing parameter interdependencies. This process attends to the natural complexity of coastal systems, particularly the random behaviour of hydrodynamic and morphological parameters [28]. The complete sampling ensures the analysis captures the full range of possible scenarios while accounting for relationships between coastal state indicators and numerical models. This considerable number of iterations enables an accurate representation of storm events during extreme conditions, particularly crucial for coastal flooding analysis. The model ensures convergence through progressive analysis of increasing sample sizes, providing confidence in the stability of the results. 33 3.1.2. Astronomical Tide Modelling The astronomical tide modelling uses data from Puertos del Estado (PORTUS) [39] oceanographic database (Table 2), specifically from the Tarragona tide gauge station. The tidal analysis revealed 22 significant constituents, implemented through the following equation (Eq. 15): 𝜂(𝑡) = 𝑍₀ + 𝛴(𝐴ᵢ 𝑐𝑜𝑠(𝜔ᵢ𝑡 + 𝜑ᵢ)) (Eq. 15) Where: - η(t): water level at time t. - Z₀: mean sea level (33.37 cm from Tarragona tide gauge)[39] . - Aᵢ: amplitude of constituent. - ωᵢ: angular frequency. - φᵢ : phase angle. Table 2: Harmonics, Tide Guge of Tarragona. Source: Puertos del Estado. https://portus.es/ 34 Tidal harmonics, including lunar (e.g., M2, N2), solar (e.g., S2), and other types (e.g., K1, O1), represent the gravitational influences of the moon, sun, and Earth’s dynamics on tidal patterns. Lunar and solar harmonics drive semi-diurnal tides with differing frequencies, while others capture diurnal cycles and complex interactions. Each harmonic is characterized by its frequency (cycles per hour), amplitude (height contribution), and phase (time offset relative to a reference point), enabling precise tide predictions. The 22 harmonic constituents are integrated into the Monte Carlo simulation to generate realistic tidal variations, which establish a continuous tidal function throughout the year. By generating random hours within the 8760 hours of a non-leap year, specific astronomical tide levels are obtained for each Monte Carlo iteration. This randomization ensures the analysis captures tidal variations across all temporal scales, from daily to seasonal cycles. Firstly, the general characteristics of the model are defined to obtain the necessary parameters for an accurate analysis. Code 1: General Characteristics As it is determined in the Code 1 text section, a random hour over a year is generated to assign determined conditions for the tide level computation. This method provides realistic astronomical tide variations for reliability analysis while maintaining the statistical independence of samples. Then, once the random hour of the year is obtained, applying the stated procedure, the following Code 2 was obtained: - Random hour generation: h ∈ [1, 8760]. - Tide level calculation: Equation 15 - Conversion to metres: results/100. %% General caracteristics total_hours = 8760; % Total hours in a non-leap year num_repetitions = 1000000; % Number of repetitions g = 9.81; % Gravitational acceleration (m/s^2) % Generation of consistent random hour selections for both analyses random_hours = randi(total_hours, num_repetitions, 1); 35 %% 2. Calculate Astronomical Tide Z0 = 33.37; % Mean sea level in cm % Harmonic components harmonics = { 'SA', 0.000114, 7.62, 249.24; 'M2', 0.080511, 3.97, 207.67; 'K1', 0.041781, 3.7, 164.87; 'O1', 0.038731, 2.4, 103.02; 'S2', 0.083333, 1.35, 229.34; 'P1', 0.041553, 1.26, 160.36; 'N2', 0.078999, 0.86, 196.71; 'S1', 0.041667, 0.66, 261.22; 'M4', 0.161023, 0.5, 346.59; 'K2', 0.083561, 0.4, 223.79; 'Q1', 0.037219, 0.31, 53.58; 'MS4', 0.163845, 0.31, 51.29; 'MN4', 0.159511, 0.19, 304.58; 'NU2', 0.079202, 0.15, 199.07; 'M3', 0.120767, 0.15, 157.55; '2N2', 0.077487, 0.14, 184.27; 'MU2', 0.077689, 0.14, 173.01; 'L2', 0.082024, 0.11, 216.42; 'T2', 0.083219, 0.1, 193.49; 'MK4', 0.164073, 0.09, 55.84; 'SK3', 0.125114, 0.08, 109.04; 'SN4', 0.162333, 0.05, 6.84; }; Code 2: Astronomical Tide Computation frequencies = cell2mat(harmonics(:, 2)); amplitudes = cell2mat(harmonics(:, 3)); phases_radians = deg2rad(cell2mat(harmonics(:, 4))); astronomical_tide = Z0 + sum(amplitudes .* cos(2 * pi * frequencies .* random_hours' + phases_radians), 1)'; astronomical_tide = astronomical_tide / 100; % Convert to metres 3.1.3. Meteorological residue The meteorological residue is a critical component in coastal flood risk analysis, particularly when assessing storm events. It accounts for water level variations caused by three key factors: 1. Atmospheric pressure changes: Low-pressure systems during storms cause water levels to rise due to the inverse barometer effect. 2. Wind stress effects: Strong winds accompanying storms can pile up water against the coastline, leading to elevated water levels and increased flood risk. 3. Large-scale oceanographic phenomena: Storms can interact with tides and ocean currents, further influencing water levels. For example, a storm coinciding with a high tide can result in considerably higher water levels. 36 Including the meteorological residual in coastal flood risk analysis ensures a more accurate assessment of potential storm surges and associated flood risks. It provides a realistic representation of the water level variations that can occur during storm events. Neglecting the meteorological residual can lead to an underestimation of flood risk, resulting in inadequate coastal protection measures, hence, incorporating this component is essential for effective coastal management to mitigate the impacts of coastal flooding. To model the meteorological residual for the Port of Tarragona, a Gaussian probability distribution is employed, with values ranging from 0.3 to 0.7. This interval is derived from empirical data presented in the Summary of Parameters Related to Sea Level and Tide that Affect Port Design and Operation Conditions specific to the port [39]. The Gaussian distribution is characterized by two key parameters: the mean (μ) and the standard deviation (σ). In this model, the mean is set to 0.5, representing the centre of the distribution, while a reduced standard deviation of 0.06 ensures that the generated values fall within the desired range. To achieve this, a truncated normal distribution approach is applied, generating random values from a standard normal distribution and transforming them to fit the specified parameters. The resulting values are filtered to retain only those within the defined bounds of 0.3 to 0.7. This approach ensures that the generated meteorological residual values are both realistic and representative of the conditions at the Port of Tarragona. These values can be utilized in coastal flood risk analysis to account for potential water level variations caused by atmospheric pressure changes, wind stress effects, and large-scale oceanographic phenomena during storm events. By fitting the Gaussian probability distribution to empirical data, the model provides a robust and reliable method for capturing the inherent variability and uncertainty of meteorological residuals, ensuring the generated values align with the specific characteristics of the Port of Tarragona. Following the outlined methodology, the Code 3 was derived: Code 3: Generation of the Meteorological Residue %% 3. Generate Meteorological Residue mu = 0.5; % Mean (center) sigma = 0.06; % Reduced standard deviation for better spread within bounds % Generate using truncated normal distribution approach z = randn(num_repetitions*2, 1); % Generate extra values for truncation z = mu + sigma * z; valid_idx = z >= 0.3 & z <= 0.7; % Find values within bounds residue_values = z(valid_idx); % Keep only valid values residue_values = residue_values(1:num_repetitions); % Take required number of samples 37 3.1.4. Wave Height Analysis and Weibull Distribution The significant wave height (Hs) at the buoy of Tarragona [39] is modelled using a Weibull probability distribution, which effectively represents both typical wave heights and extreme values. The model employs a two-regime approach to capture wave height characteristics under moderate and extreme climate conditions. Moderate Climate Regime: For the moderate climate regime, the Weibull inverse function describing the significant height of the waves in deep water is used (Eq. 16): 𝐻𝑠 = 𝐵 + 𝐴(−𝑙𝑛(1−𝑃))(1/𝐶) (Eq. 16) The parameters are as follows: - A, being the scale parameter. - B, being the location parameter. - C, being the shape parameter. Seasonal adjustments are made to the parameters to reflect varying wave conditions throughout the year: - Winter: A=0.79, B=0.22, C=1.04 - Spring: A=1.04, B=-0.06, C=1.57 - Summer: A=0.4, B=0.14, C=1.15 - Autumn: A=1.03, B=-0.01, C=1.56 The model generates random hours, the same to those used in determining astronomical tides, and random probability values between 0 and 1. These probability values are used as input to the Weibull inverse function, which assigns constants based on the respective season. It then calculates the wave height (Hs) using the Weibull inverse function. Extreme Value Regime (Hs ≥ 3m): If the computed Hs under the moderate climate regime is equal to or greater than 3 metres, which is the threshold value described in the Tarragona’s buoy for extreme events [39], the model reassigns a probability to this value and transitions to the extreme value regime. This transition involves applying adjusted constants for extreme conditions using a modified Weibull distribution. 38 Only the recalculated wave heights are considered valid, as the analysis focuses on storm events, which are the most critical and can determine the exceedance of the crest level. In this regime, a modified Weibull distribution is used to calculate the wave height (Eq. 16): 𝐻𝑠 = 𝛼+ 𝛽(−𝑙𝑛(1−𝑃))(1/𝛾) (Eq. 16) The parameters for the extreme value regime are: - α = 3.1 - β = 0.44 - γ = 0.76 The extreme value regime represents storm events, which are of particular concern for coastal and port management. It captures the wave heights that exceed the significant height limit within the average regime. This wave height analysis is based on empirical data recorded by the Tarragona Buoy, ensuring an accurate representation of the local wave climate specific to the Port of Tarragona. The primary parameters were derived from the "Average Wave Climate" dataset, while data for storm events was sourced from the document titled "Maximum Wave Extremes by Directions (Significant Wave Height)”. [39] 39 Then, because of the specified methodology, Code section 4 was developed, as follows: Code 4: Generation of Significant Wave Heights under different climate regimes. %% 4. Calculate Wave Heights and Seasonal Parameters % Define month ranges month_ranges = [ 1, 744; % January (31 days) 745, 1416; % February (28 days) 1417, 2160; % March (31 days) 2161, 2880; % April (30 days) 2881, 3624; % May (31 days) 3625, 4344; % June (30 days) 4345, 5088; % July (31 days) 5089, 5832; % August (31 days) 5833, 6552; % September (30 days) 6553, 7296; % October (31 days) 7297, 8016; % November (30 days) 8017, 8760 % December (31 days) ]; % Initialize arrays months = zeros(num_repetitions, 1); A = zeros(num_repetitions, 1); B = zeros(num_repetitions, 1); C = zeros(num_repetitions, 1); % Assign seasonal constants for i = 1:num_repetitions hour = random_hours(i); for month = 1:12 if hour >= month_ranges(month, 1) && hour <= month_ranges(month, 2) months(i) = month; % Assign constants based on season switch true case ismember(month, [12, 1, 2]) % Winter [A(i), B(i), C(i)] = deal(0.79, 0.22, 1.04); case ismember(month, 3:5) % Spring [A(i), B(i), C(i)] = deal(1.04, -0.06, 1.57); case ismember(month, 6:8) % Summer [A(i), B(i), C(i)] = deal(0.4, 0.14, 1.15); case ismember(month, 9:11) % Autumn [A(i), B(i), C(i)] = deal(1.03, -0.01, 1.56); end break; end end end % Calculate initial and adjusted Hs P = rand(num_repetitions, 1); Hs_initial = B + A .* (-log(1 - P)).^(1 ./ C); % Parameters for Hs adjustment beta = 0.44; alpha = 3.1; gamma = 0.76; adjust_indices = Hs_initial >= 3; 40 3.1.5. Directional Wave Analysis under Storm Conditions The directional wave analysis is a crucial component of the wave height model, focusing on assigning wave directions to the generated significant wave heights (Hs). The analysis employs a discrete probability distribution considering eight cardinal and intercardinal directions (N, NE, E, SE, S, SW, W, NW) to represent the primary directional patterns. The model assigns probabilities to each direction based on the observed directional distribution obtained from historical extreme events data collected by the buoy of Tarragona [39]. This data (Figure 11) provides helpful perceptions of the directional characteristics of extreme wave conditions, allowing for a more accurate representation of worst-case scenarios. To further refine the direction assignment, the model incorporates site-specific valid approach directions, accounting for local factors such as bathymetry and coastline orientation. Waves with directions outside the valid ranges are filtered out to ensure only realistic and physically possible directions are considered. It is important to note that the valid directions of incidence are defined as variables specific to each beach, and only the Hs values that follow these directions will be retained for further analysis. The assignment of wave directions to the generated Hs values, in deep water under storm events, is performed using the cumulative probability method. The probabilities for each direction are converted into cumulative probabilities, representing the likelihood of a wave coming from that direction or any preceding direction. For each filtered Hs value, a random number between 0 and 1 is generated and matched with the corresponding direction based on the cumulative probability distribution. % Initialize and calculate adjusted Hs Hs_adjusted = nan(num_repetitions, 1); P_adjusted = rand(sum(adjust_indices), 1); Hs_adjusted(adjust_indices) = alpha + beta .* (-log(1 - P_adjusted)).^(1 ./ gamma); 41 Figure 11: Wave Incidence Direction Probability. Source: Puertos del Estado, Tarragona Bouy. https://portus.es/ Hence, employing the described methodology, Code 5 this section was extracted: Code 5: Direction assignation for extreme values. %% 5. Wave Direction Analysis wave_directions = {'N', 'NE', 'E', 'SE', 'S', 'SW', 'W', 'NW'}; direction_probabilities = [4, 5, 25, 15, 24, 9, 9, 9]; % Calculate cumulative probabilities cumulative_probs = cumsum(direction_probabilities) / sum(direction_probabilities); prob_edges = [0; cumulative_probs(:)]; % Generate and assign directions P_direction = rand(sum(adjust_indices), 1); direction_indices = discretize(P_direction, prob_edges); % Convert indices to directions wave_directions_array = cell(num_repetitions, 1); wave_directions_array(adjust_indices) = wave_directions(direction_indices); % Filter valid directions wave_dirs_cat = categorical(wave_directions_array(adjust_indices)); valid_dirs_cat = categorical(valid_directions); valid_direction_mask = ismember(wave_dirs_cat, valid_dirs_cat); valid_direction_indices = false(num_repetitions, 1); valid_direction_indices(adjust_indices) = valid_direction_mask; % Create final Hs Hs_final = nan(num_repetitions, 1); Hs_final(valid_direction_indices) = Hs_adjusted(valid_direction_indices); 48 3.1.9. Convergence Analysis of Monte Carlo Simulation Results Monte Carlo simulation, while powerful for analysing coastal flooding risks, requires careful validation to ensure its results are reliable and stable. The convergence analysis serves as a critical tool for verifying that the number of simulations performed is sufficient to produce dependable probability estimates. In this study, the convergence analysis examines how the calculated probability of failure stabilizes as the number of samples increases, which means that the final probability of failure considers all the previous iterations, so the la value obtained for this probability of failure will be the most precise. The implementation follows a structured process, analysing data in increments of 50 samples per step, with a total limit of 1M iterations. To ensure the consistency of the results, a moving window of 200 samples is used to continuously evaluate stability over time [27]. The analysis employs a tolerance threshold of 0.0001 (0.01%) for determining convergence. This means that when the standard deviation of failure probabilities within the moving window falls below this threshold concerning the final probability of failure, it is considered the solution to have converged. This strict tolerance ensures high confidence in the final probability estimates. For each critical height, climate scenario (Current, RCP4.5, and RCP8.5) and year of study, the analysis tracks how the probability of failure evolves with increasing sample size. Convergence plots visualize this evolution, showing both the progression toward stable values and the point at which convergence is achieved. The convergence analysis reveals several important aspects: - First, it identifies the minimum number of simulations required to achieve reliable results, helping to optimize computational resources in future studies while ensuring accuracy. - Second, it assesses the stability of probability estimates across different scenarios. Analysing the convergence behaviour helps determine whether certain scenarios require additional simulations to achieve consistent results, indicating a higher level of variability in those conditions. 3.1.10. Reliability Analysis for Storm-Induced Coastal Flooding Coastal infrastructure faces its greatest challenges during storm events, which can generate extreme water levels and potentially catastrophic flooding. While typical sea conditions rarely threaten coastal structures, storm events can produce combinations of high waves, storm surges, and elevated water levels that pose significant risks. Therefore, reliability analysis in coastal engineering focuses primarily on these extreme conditions, as they represent the scenarios most likely to cause damage or failure. 49 The reliability index (β) quantifies a coastal system's safety level against flooding during storm events. This dimensionless parameter measures the number of standard deviations between normal operating conditions and potential failure points, expressed as (Eq. 20) [34]: 𝛽 = 𝐸[𝑀𝑆] / 𝜎[𝑀𝑆] (Eq. 20) Where: - E[MS] represents the expected value of the margin of safety. - σ[MS] denotes its standard deviation. - The margin of safety (MS) (Eq. 21) is calculated as the difference between the critical height (Hc) and the total water level (TWL) (Eq. 13) during storm conditions: 𝑀𝑆 = 𝐻𝑐 – 𝑇𝑊𝐿 (Eq. 21) For coastal flooding analysis, according to data provided in Table 3, the reliability index interpretation focuses specifically on storm event performance. Values above 3.0 indicate robust protection against severe storms, while values below 2.0 suggest vulnerability to storm-induced flooding. The minimum recommended threshold for coastal protection systems is β = 2.5, representing a balance between safety and economic feasibility [29]. Table 3: Reliability Index Range Interpretation. Source: CUR/TAW. (1990). Probabilistic design of flood defences (Report No. 141). Technical Advisory Committee on Water Defences, Center for Civil Engineering Research and Codes. Reliability Index (β) Range Interpretation for Coastal Flooding Protection During Storm Events β > 3 Robust protection against severe storms. 2.0 ≤ β ≤ 3.0 Moderate protection, meets minimum recommended threshold of β = 2.5. β < 2 Vulnerability to storm-induced flooding. Given the specific nature of storm-induced flooding, the reliability analysis employs targeted filtering to focus on potentially hazardous conditions. This filtering process considers two crucial aspects of storm events. The analysis focuses on wave directions and heights that are most relevant to storm conditions at the site. Only wave directions that align with typical storm patterns and the site's geographical exposure are considered, ensuring that the reliability assessment reflects actual storm risks. Additionally, the analysis includes only wave heights characteristic of storm events, concentrating on conditions capable of causing significant flooding. 50 By excluding minor wave states that do not pose a meaningful threat to coastal infrastructure, the reliability index provides a more accurate representation of the system's performance under truly hazardous conditions that could lead to infrastructure damage or flooding., rather than being diluted by non-threatening conditions. The methodology awards that while coastal systems may demonstrate high reliability under normal conditions, their performance during storms determines their protective value, ensuring safety remain effective against the most serious threats to coastal infrastructure. 51 Using the methodology described in the preceding three sections, the following Code 8 was developed: Code 8: Development of the sections accounting for the failure probabilities, reliability index and convergence of the obtained probabilities. %% 9. Failure, Realiability Index and Convergence with respect to Num. of Repetitions %% 9.1. Higher Critical Height %% 9.1.1. Analysis for Higher Critical Height tolerance = 0.0001; % Convergence tolerance window_size = 200; % Window size for stability step_size = 50; % Step size for sampling max_samples = num_repetitions; sample_sizes = step_size:step_size:max_samples; % Initialize arrays convergence_total_high = zeros(length(sample_sizes), length(years), length(scenarios)); convergence_points_high = zeros(length(years), length(scenarios)); failure_prob_total_high = zeros(length(years), length(scenarios)); reliability_indices_high = zeros(length(years), length(scenarios)); % Calculate base failure probability and reliability index for high critical height base_failure_high = sum(total_water_level >= critical_height_high) / num_repetitions; margin_safety_high = critical_height_high - total_water_level; beta_base_high = mean(margin_safety_high, 'omitnan') / std(margin_safety_high, 'omitnan'); % Computing Convergence, failure probabilities, and reliability indices for i = 1:length(years) for j = 1:length(scenarios) adjusted_water_level = total_water_level + slr_scenarios(i,j); % Failure probability failures = sum(adjusted_water_level >= critical_height_high); failure_prob_total_high(i,j) = failures / num_repetitions; % Reliability index margin_safety = critical_height_high - adjusted_water_level; mean_margin = mean(margin_safety, 'omitnan'); std_margin = std(margin_safety, 'omitnan'); reliability_indices_high(i,j) = mean_margin / std_margin; 52 3.1.11. Application of the model: The Ebro Delta System Figure 13: Top View of the Ebro Delta. Source: El Nacional. (2023, Date). La NASA dedica un estudio al Delta de l'Ebre y alerta de su retroceso. El Nacional. https://www.elnacional.cat/es/sociedad/nasa-dedicaestudio-delta-ebre-alerta-retrocede_688724_102.html % Convergence for k = 1:length(sample_sizes) n = sample_sizes(k); failures_n = sum(adjusted_water_level(1:n) >= critical_height_high); convergence_total_high(k,i,j) = failures_n / n; end % Find convergence point for k = window_size:length(sample_sizes) window = convergence_total_high((k-window_size+1):k,i,j); if std(window) < tolerance convergence_points_high(i,j) = sample_sizes(kwindow_size+1); break; end end end end 53 The analytical framework presented here addresses the complex coastal environment of the Ebro Delta, featuring a 50-kilometre coastline characterized by low-lying topography (Figure 13), where various natural processes interact to shape its extensive beach systems. The framework considers key environmental processes that regulate coastal behaviour throughout the delta. These include Mediterranean wave patterns and seasonal storm surge events. Sediment dynamics, historically driven by the Ebro River but now modified by upstream dam construction, play a crucial role in beach evolution. The system faces current environmental challenges, including erosion in exposed areas, sea-level rise vulnerability, and natural subsidence. This code implementation develops upon these fundamental characteristics, establishing a standardized approach for evaluating individual beaches within the delta. While maintaining consistency in the methodology, it allows for the consideration of local variations in morphodynamical conditions. The evaluation of beaches over the Ebro Delta system requires the analysis of three fundamental parameters that control their behaviour and vulnerability to coastal hazards. These parameters form the foundation of the analytical framework and are essential for understanding the dynamic response of these coastal systems. The beach slope serves as a primary indicator of coastal vulnerability and wave-beach interaction. This geometric characteristic determines how wave energy is dissipated across the beach profile and influences the potential for flooding events. For beaches in the Tarragona region and specifically the Ebro Delta, this data is obtained through the Coastal Viewer platform (Llibre Verd) from the Cartographic and Geological Institute of Catalonia [41]. By selecting the specific beach of interest, detailed slope information can be extracted, providing a complete understanding of the beach profile configuration. Wave incidence directions constitute the second critical parameter in the analysis framework. These directional components describe the predominant angles at which waves approach the coastline. The analysis of wave directions is conducted using cartographic tools, where a reference line is established for the beach of interest. From the eight possible cardinal and intercardinal directions, five primary wave approach directions are identified as potentially affecting the beach. This directional distribution significantly influences sediment transport patterns along the shore and determines the probability of wave impact directions on different beaches. The third essential parameter is the critical wave height, which represents a crucial threshold parameter in beach flooding analysis, particularly significant in low-lying areas such as the Ebro Delta. This parameter, obtained through the VISSIR [40] visualization platform of the Cartographic and Geological Institute of Catalonia, helps identify the highest points that serve as the last line of defence against flooding. When floodwater breaches this highest point, it can rapidly flow to lower elevations and inland areas, potentially causing beach erosion and infrastructure damage. 54 However, a comprehensive flooding assessment requires consideration of both the highest elevations and lower-lying areas. While the highest point acts as a critical threshold for flooding protection, lower-lying sections can function as channels directing water toward inner areas, even before the highest point is overtopped. This combined analysis of both high and low points, all accessible through VISSIR, enables a detailed understanding of flooding mechanisms and vulnerabilities. Together, these three parameters establish a comprehensive framework for analysing beach behaviour within the Ebro Delta system. The Marquesa and The Trabucador beaches were specifically selected for this analysis due to their contrasting positions within the delta, each located on opposite sides of the system, demonstrating how local characteristics influence beach behaviour and vulnerability to coastal processes. The Marquesa, located in the northern hemidelta, and The Trabucador, positioned in the southern section, experience distinct wave patterns and morphodynamical conditions. 3.1.12. Beach Characteristics: Site-Specific Analysis for The Marquesa Beach Located in the northern hemidelta of the Ebro Delta system (Figure 14), Marquesa Beach (Figure 15) represents a coastal environment with a characteristic profile of fine golden sand composition, creating a moderate slope that influences wave run-up processes. Its position makes it particularly susceptible to northerly wind patterns, specifically the Tramontana, which plays a significant role in sediment transport dynamics. The beach system incorporates a well-developed berm and experiences relatively shallow waters, contributing to its specific wave transformation patterns. The backbeach area includes both natural and modified zones, measuring critical wave heights is particularly relevant for flood risk management. Its orientation allows for wave approach from multiple directions, though it remains partially sheltered from certain wave angles by its position within the delta configuration, however in this study, as the input data was obtained from buoys located in deep water all possible directions were considered, to consider the all the actual probable cases. 55 Figure 14: Expanded top view of The Marquesa. Source: Google maps. Figure 15: Top View and Orientation of the Marquesa Beach: Source: Google Maps. 56 - The Marquesa Beach Slope Analysis The Marquesa beach slope value is obtained through the Coastal Viewer platform (Llibre Verd) from the Cartographic and Geological Institute of Catalonia [41]. The platform provides detailed topographic and bathymetric data that allows for precise slope calculations. Analysis of the beach profile shows that The Marquesa has a slope of 0.05 (5%). This measurement reflects the beach's characteristic gentle gradient, which is typical of its position in the northern hemidelta. The slope calculation considers the profile from the shoreline extending offshore, providing an average gradient that is representative of the beach face used in wave run-up calculations. - The Marquesa Beach Wave Exposure Given The Marquesa beach's location in the Ebro Delta, as it is determined in Figure 15, waves can approach from the following principal directions: NW (Northwest), N (North), NE (Northeast), E (East), and SE (Southeast). The most significant wave action typically comes from the easterly (E) and northeasterly (NE) directions, as these waves arrive directly from the Mediterranean Sea. The beach also experiences wave patterns from the southeastern (SE) direction, particularly during certain meteorological conditions. - Critical Height Analysis for The Marquesa Beach For The Marquesa Beach, the critical height has been established at 2.2 metres, representing the highest elevation point that serves as the last defence against flooding. This threshold was determined based on the beach's topographic characteristics, including its moderate slope and relatively wide beach profile, which provides some natural protection against wave action. However, a comprehensive flooding analysis reveals that the inner areas become vulnerable at a much lower threshold of 1.3 metres, where lower-lying sections can function as channels directing water toward agricultural areas and infrastructure. Through topographic data obtained from VISSIR (Figure 16), it is possible to observe how these different elevation points interact during flooding events. When waves exceed the lower 1.3-metre threshold, water can begin infiltrating through these natural channels, potentially affecting inland areas even before the highest point is overtopped. Once waves surpass the 2.2-metre critical height, the entire beach system becomes vulnerable to more severe flooding, with rapid inundation of lower areas and potential damage to nearby infrastructure. 57 Figure 16: Critical height and low elevations for the Marquesa Beach . Source: Institut Cartogràfic i Geològic de Catalunya. VISSIR3 - Visualizador de información geográfica. http://srv.icgc.cat/vissir3/ This dual threshold consideration demonstrates the importance of analysing both high and low points in flood risk assessment. While the 2.2-metre critical height provides a conservative measure for maximum protection, the lower 1.3-metre threshold alerts us to initial flooding risks through natural channels. This comprehensive understanding enables more effective flood risk management strategies, which are particularly crucial for protecting the adjacent agricultural areas and infrastructure that characterize Marquesa Beach's back-beach environment. In this way, the local characteristics of the Marquesa are defined in Code 9. Code 9: Definition of the morphological characteristics and wave exposure directions for the Marquesa. %% The Marquesa slope = 0.05; % Beach slope critical_height_high = 2.2; % Higher critical water level (m) critical_height_low = 1.3; % Lower critical water level (m) valid_directions = {'NW', 'N', 'NE', 'E', 'SE'}; 64 4. Results The results section presents the analytical outcomes obtained through the previously described modelling framework. The analysis focuses on two distinct coastal locations: The Marquesa and The Trabucador beaches, where multiple visualization sets have been generated to illustrate various aspects of coastal vulnerability under different scenarios. The assessment involves a complete analysis of site-specific responses generated by the Monte Carlo simulation model. The model was executed for 1 million iterations using MATLAB software on an HP Pavilion computer, with each site requiring nearly one hour of computation. While maintaining consistent computational parameters and modifying the characteristics of each location. The analysis framework enables direct attribution of observed variations to the distinctive morphological features and directional exposure patterns unique to each beach location. The analytical results demonstrate water level projections under various climate change scenarios, incorporating both near-term and long-term forecasts, while accounting for local subsidence effects. These projections are complemented by detailed probability of failure assessments, which have been evaluated through two distinct approaches: considering the complete dataset of samples and specifically focusing on storm event conditions. Additionally, the analysis includes convergence studies and a comprehensive reliability assessment that evaluates safety levels against flooding during storm events, providing quantitative measures of each beach's resilience under extreme conditions. 65 4.1. The Marquesa 4.1.1. Wave analysis Figure 21: Hs Distribution for the Marquesa in function of the climate. Source: Main. The figure shown above (Figure 21) presents a comparative analysis of wave height distributions at The Marquesa beach, illustrating the impact of directional and storm events filtering criteria on significant wave heights (Hs). The Initial Hs Distribution contains the complete wave height dataset, ranging from 0 to 10 metres, with most of the wave heights concentrated between 0 and 2 metres. This distribution represents all recorded wave conditions, achieving 100% valid samples, and demonstrates more values on the left and a long tail on the right, which is a pattern typical of wave height distributions in the Mediterranean region. The Final Hs Distribution focuses exclusively on storm events, representing only 0.5% of the total recorded cases. These events are specifically filtered to include waves approaching from NW, N, NE, E, and SE directions with heights equal to or exceeding 3 metres from the Initial Hs Distribution. This filtered dataset has been recalculated to analyse extreme conditions within these specific directional constraints, with wave heights primarily ranging between 3 and 6 metres. The resulting distribution is afterwards utilized for run-up computations, ensuring that the analysis captures the most critical conditions that influence coastal flooding risks at The Marquesa. 66 4.1.2. Water Level Components Figure 22: Water Level Components for The Marquesa beach. Source: Main. The figure showed above (Figure 22), illustrates the key components contributing to water level variations in current scenario at The Marquesa beach through four distinct distributions. The Astronomical Tide Distribution shows a pattern ranging from 0.17 to 0.53 metres, with peak values concentrated around 0.34 metres, related with the mean sea level of 33.37cm. The Meteorological Residue Distribution describes a Gaussian function spanning from 0.30 to 0.70 metres, centred at 0.50 metres, reflecting the atmospheric pressure and wind influences on water levels. Wave Run-up (R2%), considering Hs filtered for storm events values in function of their direction, exhibits again a right-skewed distribution with most events concentrated between 0.5 and 0.8 metres, declining toward higher values. The Total Water Level Distribution, integrating all components, ranges from 1.2 to 2.6 metres. Notably, a significant portion of the distribution exceeds the Low Critical Height (1.30m), while values above the High Critical Height (2.20m) are infrequent, suggesting moderate vulnerability to extreme water level events at The Marquesa beach. 67 4.1.3. Failure Probability during Storm Events Figure 23: Failure Probabilities In Case of Storm Event by Critical Height (the Marquesa). Source: Main. Figure 23, illustrates the failure probabilities (0 to 1) for The Marquesa beach during storm events across different climate change scenarios and time horizons, analysing both the high (2.2m) and low (1.3m) critical heights. For the high critical height, representing the beach's elevated areas, the current scenario shows minimal failure probability (0.002) in 2025. Of course, probability values appear identical for all three scenarios (Current, RCP4.5, and RCP8.5) in 2025, as the effects of different climate change trajectories have not yet had time to materialize in this nearterm timeframe. This risk remains relatively low in 2050, with slight increases to 0.022 and 0.030 for RCP4.5 and RCP8.5 scenarios, respectively. However, by 2100, the probability of failure increases substantially, reaching 0.440 under RCP4.5 and 0.999 under RCP8.5, indicating a significant long-term vulnerability when accounting for subsidence effects and climate change scenarios. Regarding the low critical height, corresponding to the beach's lower sections, the analysis reveals consistently high failure probabilities across all scenarios and timeframes. Starting with 0.932 in 2025 for all scenarios, the probability increases to 0.991-1.000 by 2050 and maintains these elevated levels through 2100. These results demonstrate that during storm events, the lower sections of The Marquesa Beach face immediate and persistent flooding risks, with climate change and subsidence effects amplifying this vulnerability over time. 68 The explicit contrast between high and low critical heights emphasizes the spatial variability of flooding risks within Marquesa Beach, while the temporal progression across scenarios highlights the increasing impact of climate change and subsidence on coastal flooding probabilities. 4.1.4. Water Level Exceedance The water level exceedance probability graph for Marquesa Beach illustrates an increasing flood risk pattern over time (Figure 24). Current conditions in 2025 show relatively low probabilities of exceeding the high critical height (2.2m). However, both moderate (RCP4.5) and severe (RCP8.5) climate change scenarios indicate progressively higher water levels by 2050, with substantial increases projected for 2100. The low critical height (1.3m) remains consistently vulnerable across all timeframes and scenarios. Most notably, the RCP8.5 scenario for 2100 demonstrates a significant probability of exceedingly even the high critical threshold, indicating severe future flooding risks without intervention. Figure 24: Water Level Exceedance In Case of Storm Event by Critical Height (the Marquesa). Source: Main. 4.1.5. Convergence Analysis This section develops the convergence analysis of the failure probability with respect to all the samples studied. The Marquesa beach's high critical height (2.20m) reveals distinct patterns across different climate scenarios, with convergence determined through specific criteria: a tolerance of 0.0001, a stability window of 200 samples, and progressive evaluation at 50-sample intervals (Figure 24). 69 Figure 25: Convergence Analysis-Current Scenario, High Critical Height (the Marquesa). Source: Main. Under the current scenario, as shown on the figure above (Figure 25), the analysis demonstrates efficient convergence with minimal failure probability. The model achieves stability at n=50 samples, meeting the convergence criteria with consistent evaluations across the specified window size, while maintaining a notably low failure probability (Pf=0.0000) through 2025 and 2050. For 2100, a slight increase in failure probability is observed (Pf=0.0001), though still converging efficiently within the same sampling parameters. Figure 26: Convergence Analysis-RCP4.5 Scenario, High Critical Height (the Marquesa). Source: Main. 70 The RCP4.5 scenario, shown on Figure 26, exhibits more complex convergence behaviour. While 2025 and 2050 maintain quick convergence at n=50 samples with low failure probabilities (Pf=0.0000 and 0.0001 respectively), the 2100 projection requires substantially more samples (n=11700) to satisfy the convergence criteria. This increased sample requirement, yielding a failure probability of 0.0022, reflects the growing uncertainty in long-term projections, necessitating more iterations to achieve stability within the specified tolerance. Figure 27: Convergence Analysis-RCP8.5 Scenario, High Critical Height (the Marquesa). Source: Main. The RCP8.5 scenario (Figure 27) presents the most challenging convergence patterns. While early periods maintain efficient convergence like other scenarios, the 2100 projection demands the highest number of samples (n=11900) to reach stability within the defined tolerance and window size, resulting in a failure probability of 0.0049. The convergence pattern shows initial fluctuations before stabilizing across the specified window, indicating increased complexity in modelling more severe climate change impacts under the strict convergence criteria established. This progressive increase in required samples and failure probabilities across scenarios and timeframes demonstrates the growing complexity and uncertainty in modelling more severe climate change impacts, particularly for long-term projections, while ensuring robust results through rigorous convergence criteria. The convergence analysis for Marquesa Beach's low critical height (1.30m) demonstrates a notably different pattern compared to the high critical height, with consistent behaviour across all climate scenarios while adhering to the established convergence criteria (tolerance of 0.0001, stability window of 200 samples, and 50sample intervals). 71 Figure 28: Convergence Analysis-Current Scenario, Low Critical Height (the Marquesa). Source: Main. Figure 29: Convergence Analysis-RCP4.5 Scenario, Low Critical Height (the Marquesa). Source: Main. 72 Figure 30: Convergence Analysis-RCP8.5 Scenario, Low Critical Height (the Marquesa). Source: Main. For all scenarios (Current, RCP4.5, and RCP8.5), the analysis requires a significantly higher number of samples (n=11900) to achieve convergence, indicating greater complexity in modelling the lower elevation behaviour. The failure probabilities show a progressive pattern, starting at Pf=0.0046 for 2025 and increasing slightly to Pf=0.0049 for both 2050 and 2100 (Figure 28, 29 and 30, respectively). The convergence patterns exhibit initial fluctuations before stabilizing, particularly evident in the first 100000 samples. After this point, all scenarios demonstrate remarkable stability, maintaining consistent failure probabilities across the remainder of the simulation period. This stability is particularly notable as it persists across different climate scenarios, suggesting that the low-lying areas' vulnerability is more influenced by their topographic position than by the specific climate change scenario. The succeeding convergence analysis examines flooding under critical storm conditions, focusing on waves that exceed height thresholds and approach from key directions. This assessment of extreme events demonstrates the Monte Carlo simulation's robustness in evaluating beach performance during severe conditions. The model achieves stable probability calculations even when analysing this specific subset of hazardous scenarios, validating its effectiveness for coastal vulnerability assessment. While the higher elevations showed increasing complexity and higher failure probabilities with more severe climate scenarios, the lower elevations maintain consistent behaviour, indicating their persistent vulnerability regardless of the climate change trajectory, suggesting that low-lying areas require immediate attention regardless of the climate change scenario considered (Figures 31 and 32). 73 Figure 31: Convergence Analysis-Valid Samples High Critical Height (the Marquesa). Source: Main Figure 32: Convergence Analysis-Valid Samples Low Critical Height (the Marquesa). Source: Main. 80 4.2.4. Water Level Exceedance This water level exceedance probability graph for Trabucador Beach reveals concerning vulnerability patterns. Under current conditions in 2025, the beach already shows a notable probability of exceeding its high critical height (2.1m). The progression through time demonstrates a marked increase in flooding risk, with the 2050 projections showing significant rightward shifts in the probability curves for both moderate (RCP4.5) and severe (RCP8.5) climate change scenarios (Figure 38). The situation becomes particularly severe by 2100, where the probability curves indicate frequent exceedance of critical thresholds. Most notably, the RCP8.5 scenario for 2100 shows extremely high probabilities of water levels exceeding both critical heights, suggesting near-certain flooding conditions during severe weather events. The low critical height (1.2m) appears consistently vulnerable across all scenarios, with 100% exceedance probability indicated by the flat upper portions of all curves. Comparing this to the Marquesa beach, the Trabucador demonstrates greater overall vulnerability, with higher exceedance probabilities across all scenarios and periods. This increased vulnerability shows clearly how all the probability curves have shifted rightward, indicating that higher water levels occur more frequently at this location. The analysis suggests that the Trabucador requires more immediate attention for coastal protection measures, particularly given its rapid deterioration in protective capacity under future climate scenarios. Figure 38: Water Level Exceedance In Case of Storm Event by Critical Height (the Trabucador). Source: Main. 81 4.2.5. Convergence Analysis The convergence analysis for Trabucador Beach's high critical height (2.10m) reveals distinct patterns across climate scenarios and temporal horizons while following the established convergence criteria (tolerance of 0.0001, stability window of 200 samples, and 50-sample steps). Figure 39: Convergence Analysis-Current Scenario, High Critical Height (the Trabucador). Source: Main. Under the current scenario (Figure 39), the analysis demonstrates a progressive increase in required samples for convergence. For 2025, convergence is achieved efficiently at n=50 samples with a failure probability of 0.0002. The 2050 projection requires n=1800 samples to reach stability with Pf=0.0003, while 2100 demands n=8650 samples to achieve convergence with Pf=0.0007, indicating increasing complexity in long-term projections. 82 Figure 40: Convergence Analysis-RCP4.5 Scenario, High Critical Height (the Trabucador). Source: Main. The RCP4.5 scenario (Figure 40) shows more pronounced changes in convergence behaviour. While maintaining efficient convergence for 2025 (n=50, Pf=0.0002), the 2050 projection requires n=5100 samples to achieve stability with Pf=0.0013. The 2100 horizon exhibits the highest complexity, needing n=16650 samples to converge with Pf=0.0066, reflecting increased uncertainty in long-term climate change impacts. Figure 41: Convergence Analysis-RCP8.5 Scenario, High Critical Height (the Trabucador). Source: Main. 83 The RCP8.5 scenario shown in Figure 41, maintains the same convergence characteristics for 2025 and 2100 as RCP4.5, but shows distinct behaviour for 2050, requiring n=9050 samples for convergence with Pf=0.0018. This increased sample requirement for mid-century projections suggests complexity in modelling the transition period under more severe climate change conditions. The convergence analysis for Trabucador Beach's low critical height (1.20m) demonstrates consistent behaviour across all climate scenarios and temporal horizons, following the established convergence criteria (tolerance of 0.0001, stability window of 200 samples, and 50-sample steps). Figure 42: Convergence Analysis-Current Scenario, Low Critical Height (the Trabucador). Source: Main. Figure 43: Convergence Analysis-RCP4.5 Scenario, Low Critical Height (the Trabucador). Source: Main. 84 Figure 44: Convergence Analysis-RCP8.5 Scenario, Low Critical Height (the Trabucador). Source: Main. Notably, all scenarios (Current, RCP4.5, and RCP8.5) and timeframes (2025, 2050, and 2100) exhibit identical convergence characteristics, requiring n=16650 samples to achieve stability with a failure probability of 0.0066. This consistency suggests that the low-lying areas of The Trabucador maintain similar vulnerability patterns regardless of the climate change scenario or temporal horizon considered (Figures 42, 43 and 44, respectively). The convergence patterns show initial fluctuations in the first 100000 samples before stabilizing, with all scenarios maintaining steady failure probabilities the simulation period. This stability is remarkable as it persists across different climate scenarios and timeframes. By analysing exclusively storm event conditions, to determine the probability of failure, where wave heights exceed established thresholds and approach from valid directions, the convergence analysis still demonstrates robust statistical stability. The analysis confirms that even when constraining the dataset to these extreme events, the Monte Carlo simulation achieves reliable convergence, validating the methodology's effectiveness in evaluating beach performance during storms. So, the convergence assessment for this subset of data shows that the probability of failure calculations maintains stability despite the reduced sample size. 85 Figure 45: Convergence Analysis-Valid Samples High Critical Height (the Trabucador). Source: Main Figure 46: Convergence Analysis-Valid Samples Low Critical Height (the Trabucador). Source: Main 86 Similar to what happened in The Marquesa, higher elevations show increasing failure probabilities under more severe climate scenarios, which stabilize after the first 100000 samples, as shown in Figure 45, while The Trabucador's low-lying areas demonstrate the worst uniform behaviour (Figure 46). This suggests that The Trabucador's lower elevations present more consistent vulnerability characteristics, maintaining steady failure probabilities across all scenarios and timeframes. This contrasts with The Marquesa's more variable convergence patterns and generally lower failure probabilities, particularly in earlier timeframes and less severe climate scenarios. 4.2.6 Reliability Analysis Figure 47: Reliability Indices by Scenario. High Critical Height for the Trabucador. Source: Main. The reliability indices for Trabucador Beach reveal awareness of its vulnerability to flooding, with distinct patterns emerging across different elevations and time horizons (Figure 47). For the high critical height (2.10m), the initial moderate reliability index of 2.45 in 2025 indicates that even the elevated areas of The Trabucador start with only acceptable safety levels against flooding during storm events. This contrasts sharply with The Marquesa's robust initial index of 4.83, suggesting that The Trabucador's higher elevations are inherently more vulnerable. 87 The rapid deterioration by 2050, particularly under climate change scenarios (dropping to 0.59 and 0.31), indicates a transition toward unsafe conditions much earlier than The Marquesa. By 2100, the negative indices (-1.70 and -3.41 under RCP scenarios) signify that even these higher elevations will likely experience regular flooding during storms, representing a fundamental shift in the beach's protective capacity. Figure 48: Reliability Indices by Scenario. Low Critical Height for the Trabucador. Source: Main. The low critical height (1.20m) analysis reveals even more concerning implications. The strongly negative initial indices around -2.70 in 2025 indicate that these areas are already in a critical state, with flooding being a near-certainty during storm events. The progression to extreme negative values by 2100 (-6.85 and -8.56 under RCP scenarios) suggests these lower areas will become essentially uninhabitable during storm conditions, requiring immediate and significant adaptation strategies (Figure 48). The contrast between The Trabucador and The Marquesa's reliability indices emphasizes how local coastal characteristics significantly influence vulnerability. The Trabucador's more severe indices across all scenarios can be attributed to its specific morphological features, including its wider directional exposure to waves and steeper beach slope. These characteristics make it not only more susceptible to immediate flooding risks but also more sensitive to the combining effects of sea-level rise and climate change, suggesting a need for more urgent and comprehensive adaptation strategies compared to The Marquesa. 88 5. Conclusions Based on the comprehensive analysis of The Marquesa and The Trabucador beaches in the Ebro Delta, several key insights emerge concerning their effectiveness as natural coastal defences. The findings reveal distinct patterns in how beach morphology and orientation influence protective capacity, with significant implications for coastal management strategies. The research establishes that geometric characteristics fundamentally determine a beach's defensive capabilities. While both study sites experience identical tidal and meteorological conditions, their responses to wave action differ significantly due to their distinct configurations. The Marquesa's higher elevations currently demonstrate robust protection against flooding, though this effectiveness is projected to diminish gradually through mid-century before facing substantial challenges by 2100 under severe climate scenarios. In contrast, The Trabucador exhibits increased sensitivity to environmental forces even at its elevated sections, with protection levels starting at acceptable but concerning values and confirming accelerated deterioration under climate change projections. The analysis identifies beach elevation as the critical determinant of flood protection capability. The lower sections of both beaches already demonstrate significant vulnerability and show persistent susceptibility to flooding across all climate scenarios, indicating that elevation challenges constitute a fundamental concern independent of climate change trajectories. Even though, The Trabucador's situation is particularly even more serious than the Marquesa’s, which requires more immediate attention for coastal protection measures, particularly given its rapid deterioration in protective capacity under future climate scenarios. The methodology's reliability is validated through extensive statistical testing, demonstrating consistent results across varying sample sizes while accounting for the increasing complexity of long-term climate projections. This robust analytical framework provides faithful evaluations for both current and future scenarios, though its validity depends on maintaining assumed beach geometries through regular maintenance programs. When extended to the general Ebro Delta system, the findings suggest that beaches with similar physical characteristics will encounter comparable challenges. Northern-facing beaches benefit from natural protection against predominant storm directions, while those with broader directional exposure face more immediate risks requiring frequent intervention. The consistent vulnerability patterns observed at lower elevations throughout the study area emphasize the need for immediate attention to low-lying sections across the delta. 89 These conclusions underscore the urgent requirement for a comprehensive coastal management strategy that integrates both maintenance and adaptation components. Areas below 1.30 metres demonstrate immediate vulnerability regardless of location, while the accelerating deterioration of protective capacity under climate change scenarios demands adaptive management considerations. This strategy must address both immediate vulnerabilities and long-term climate impacts through continuously adjusted maintenance programs. Looking ahead, these findings establish a scientific foundation for evidence-based coastal protection strategies. The analysis demonstrates that maintaining current beach geometries alone will not be enough, strategic interventions to enhance protective features, particularly through elevation modifications, will be essential for long-term coastal resilience. As climate change continues to reshape coastal dynamics, the implementation of these insights becomes increasingly critical for preserving the Ebro Delta's natural and cultural heritage. The message from coastal regions grows clearer and more urgent each day. As waters rise and storms intensify, the early signs of climate change's impact on coastlines become increasingly apparent. The Ebro Delta's beaches tell a story that coastal communities worldwide need to hear, they serve as an early warning system, demonstrating the consequences when natural defences begin to fail. The research leaves no room for doubt, without immediate action to protect and strengthen these natural barriers, the losses extend far beyond sand and shoreline. At concern are established communities, rich ecosystems that have been there for centuries, and cultural heritage that, once claimed by the sea, remains irrecoverable. This moment stands as a critical crisis. The advancing waves will not pause for extended deliberation or delayed decisions. The scientific evidence stands clear, the data is in front of us, and the path forward demands immediate action. Tomorrow's rising tides advance inevitably, and the resilience of coastal regions depends entirely on actions taken in the immediate term. 96 [44] United Nations. (2015). Goal 13: Take urgent action to combat climate change and its impacts. In Transforming our world: The 2030 agenda for sustainable development. Retrieved from https://sdgs.un.org/goals/goal13 [45] International Energy Agency. (2019). Global energy & CO2 status report 2019. https://www.iea.org/reports/global-energy-co2-status-report-2019 97 Appendices 1. Model code clear all; close all; %% 1. Characteristics of the model %load('base_monte_carlo.mat'); load('my_trabucador.mat') %% The Marquesa % slope = 0.05; % Beach slope % critical_height_high = 2.2; % Higher critical water level (m) % critical_height_low = 1.3; % Lower critical water level (m) % valid_directions = {'NW', 'N', 'NE', 'E', 'SE'}; %% The Trabucador slope = 0.065; % Beach slope critical_height_high = 2.1; % Higher critical water level (m) critical_height_low = 1.2; % Lower critical water level (m) valid_directions = {'NE', 'E', 'SE','S','SW'}; %% General caracteristics total_hours = 8760; % Total hours in a non-leap year num_repetitions = 1000000; % Number of repetitions g = 9.81; % Gravitational acceleration (m/s^2) % Generate consistent random hour selections for both analyses random_hours = randi(total_hours, num_repetitions, 1); %% 2. Calculate Astronomical Tide Z0 = 33.37; % Mean sea level in cm % Harmonic components harmonics = { 'SA', 0.000114, 7.62, 249.24; 'M2', 0.080511, 3.97, 207.67; 'K1', 0.041781, 3.7, 164.87; 'O1', 0.038731, 2.4, 103.02; 'S2', 0.083333, 1.35, 229.34; 'P1', 0.041553, 1.26, 160.36; 'N2', 0.078999, 0.86, 196.71; 'S1', 0.041667, 0.66, 261.22; 'M4', 0.161023, 0.5, 346.59; 'K2', 0.083561, 0.4, 223.79; 'Q1', 0.037219, 0.31, 53.58; 'MS4', 0.163845, 0.31, 51.29; 'MN4', 0.159511, 0.19, 304.58; 'NU2', 0.079202, 0.15, 199.07; 'M3', 0.120767, 0.15, 157.55; '2N2', 0.077487, 0.14, 184.27; 'MU2', 0.077689, 0.14, 173.01; 'L2', 0.082024, 0.11, 216.42; 'T2', 0.083219, 0.1, 193.49; 'MK4', 0.164073, 0.09, 55.84; 'SK3', 0.125114, 0.08, 109.04; 'SN4', 0.162333, 0.05, 6.84; }; 98 [A(i), B(i), C(i)] = deal(1.03, -0.01, 1.56); frequencies = cell2mat(harmonics(:, 2)); amplitudes = cell2mat(harmonics(:, 3)); phases_radians = deg2rad(cell2mat(harmonics(:, 4))); astronomical_tide = Z0 + sum(amplitudes .* cos(2 * pi * frequencies .* random_hours' + phases_radians), 1)'; astronomical_tide = astronomical_tide / 100; % Convert to metres %% 3. Generate Meteorological Residue mu = 0.5; % Mean (center) sigma = 0.06; % Reduced standard deviation for better spread within bounds % Generate using truncated normal distribution approach z = randn(num_repetitions*2, 1); % Generate extra values for truncation z = mu + sigma * z; valid_idx = z >= 0.3 & z <= 0.7; % Find values within bounds residue_values = z(valid_idx); % Keep only valid values residue_values = residue_values(1:num_repetitions); % Take required number of samples %% 4. Calculate Wave Heights and Seasonal Parameters % Define month ranges month_ranges = [ 1, 744; % January (31 days) 745, 1416; % February (28 days) 1417, 2160; % March (31 days) 2161, 2880; % April (30 days) 2881, 3624; % May (31 days) 3625, 4344; % June (30 days) 4345, 5088; % July (31 days) 5089, 5832; % August (31 days) 5833, 6552; % September (30 days) 6553, 7296; % October (31 days) 7297, 8016; % November (30 days) 8017, 8760 % December (31 days) ]; % Initialize arrays months = zeros(num_repetitions, 1); A = zeros(num_repetitions, 1); B = zeros(num_repetitions, 1); C = zeros(num_repetitions, 1); % Assign seasonal constants for i = 1:num_repetitions hour = random_hours(i); for month = 1:12 if hour >= month_ranges(month, 1) && hour <= month_ranges(month, 2) months(i) = month; % Assign constants based on season switch true case ismember(month, [12, 1, 2]) % Winter [A(i), B(i), C(i)] = deal(0.79, 0.22, 1.04); case ismember(month, 3:5) % Spring [A(i), B(i), C(i)] = deal(1.04, -0.06, 1.57); case ismember(month, 6:8) % Summer [A(i), B(i), C(i)] = deal(0.4, 0.14, 1.15); case ismember(month, 9:11) % Autumn 99 Tp(valid_indices) = 4.32 * Hs_final(valid_indices).^0.43; end break; end end end % Calculate initial and adjusted Hs P = rand(num_repetitions, 1); Hs_initial = B + A .* (-log(1 - P)).^(1 ./ C); % Parameters for Hs adjustment beta = 0.44; alpha = 3.1; gamma = 0.76; adjust_indices = Hs_initial >= 3; % Initialize and calculate adjusted Hs Hs_adjusted = nan(num_repetitions, 1); P_adjusted = rand(sum(adjust_indices), 1); Hs_adjusted(adjust_indices) = alpha + beta .* (-log(1 - P_adjusted)).^(1 ./ gamma); % Count the number of valid (non-NaN) Hs_adjusted values num_valid_Hs_adjusted = sum(~isnan(Hs_adjusted)); percentile_temporal=(num_valid_Hs_adjusted/num_repetitions)*100; %% 5. Wave Direction Analysis wave_directions = {'N', 'NE', 'E', 'SE', 'S', 'SW', 'W', 'NW'}; direction_probabilities = [4, 5, 25, 15, 24, 9, 9, 9]; % Calculate cumulative probabilities cumulative_probs = cumsum(direction_probabilities) / sum(direction_probabilities); prob_edges = [0; cumulative_probs(:)]; % Generate and assign directions P_direction = rand(sum(adjust_indices), 1); direction_indices = discretize(P_direction, prob_edges); % Convert indices to directions wave_directions_array = cell(num_repetitions, 1); wave_directions_array(adjust_indices) = wave_directions(direction_indices); % Filter valid directions wave_dirs_cat = categorical(wave_directions_array(adjust_indices)); valid_dirs_cat = categorical(valid_directions); valid_direction_mask = ismember(wave_dirs_cat, valid_dirs_cat); valid_direction_indices = false(num_repetitions, 1); valid_direction_indices(adjust_indices) = valid_direction_mask; % Create final Hs Hs_final = nan(num_repetitions, 1); Hs_final(valid_direction_indices) = Hs_adjusted(valid_direction_indices); %% 6. Calculate Wave Parameters % Calculate Tp Tp = nan(num_repetitions, 1); valid_indices = valid_direction_indices & ~isnan(Hs_final); 100 % Calculate deep water wavelength (L0) L0 = nan(num_repetitions, 1); L0(valid_indices) = (g * Tp(valid_indices).^2) ./ (2 * pi); % Calculate Irribarren number (xi) xi = nan(num_repetitions, 1); xi(valid_indices) = slope ./ sqrt(Hs_final(valid_indices) ./ L0(valid_indices)); %% 7. Calculate Wave Run-up (R2%) R2 = nan(num_repetitions, 1); % For xi < 3 low_xi_indices = valid_indices & (xi < 3); R2(low_xi_indices) = 0.73 * slope * sqrt(Hs_final(low_xi_indices) .* L0(low_xi_indices)); % For xi >= 3 high_xi_indices = valid_indices & (xi >= 3); R2(high_xi_indices) = 1.1 * (0.35 * slope * sqrt(Hs_final(high_xi_indices) .* L0(high_xi_indices)) + ... sqrt(Hs_final(high_xi_indices) .* L0(high_xi_indices) .* (0.563 * slope^2 + 0.004)) / 2); total_water_level = astronomical_tide + residue_values + R2; %% 8. Sea Level Rise (SLR) + Subcidence Scenario % Sea Level Rise Scenarios (in metres) % Moderate scenario (RCP4.5) slr_2050_moderate = 0.25; % 25 cm rise by 2050 slr_2100_moderate = 0.50; % 50 cm rise by 2100 % High emission scenario (RCP8.5) slr_2050_high = 0.30; % 30 cm rise by 2050 slr_2100_high = 0.80; % 80 cm rise by 2100 % Local subsidence rate (metres/year) subsidence_rate = 0.003; % 3 mm/year % Time periods for analysis years = [2025, 2050, 2100]; scenarios = {'Current', 'RCP4.5', 'RCP8.5'}; % Calculate SLR for different scenarios and years slr_scenarios = zeros(length(years), length(scenarios)); for i = 1:length(years) year = years(i); years_from_now = year - 2025; % Current scenario (only subsidence) slr_scenarios(i, 1) = subsidence_rate * years_from_now; % RCP4.5 (moderate) scenario if year <= 2050 slr_scenarios(i, 2) = (slr_2050_moderate * (year - 2025)/(2050 - 2025)) + ... (subsidence_rate * years_from_now); 101 else slr_scenarios(i, 2) = slr_2050_moderate + ... (slr_2100_moderate - slr_2050_moderate) * (year - 2050)/(2100 - 2050) + ... (subsidence_rate * years_from_now); end % RCP8.5 (high) scenario if year <= 2050 slr_scenarios(i, 3) = (slr_2050_high * (year - 2025)/(2050 - 2025)) + ... (subsidence_rate * years_from_now); else slr_scenarios(i, 3) = slr_2050_high + ... (slr_2100_high - slr_2050_high) * (year - 2050)/(2100 - 2050) + ... (subsidence_rate * years_from_now); end end %% 9. Failure, Realiability Index and Convergence with respect to Num. of Repetitions %% 9.1. Higher Critical height %% 9.1.1. Analysis for Higher Critical Height tolerance = 0.0001; % Convergence tolerance window_size = 200; % Window size for stability step_size = 50; % Step size for sampling max_samples = num_repetitions; sample_sizes = step_size:step_size:max_samples; % Initialize arrays convergence_total_high = zeros(length(sample_sizes), length(years), length(scenarios)); convergence_points_high = zeros(length(years), length(scenarios)); failure_prob_total_high = zeros(length(years), length(scenarios)); reliability_indices_high = zeros(length(years), length(scenarios)); % Calculate base failure probability and reliability index for high critical height base_failure_high = sum(total_water_level >= critical_height_high) / num_repetitions; margin_safety_high = critical_height_high - total_water_level; beta_base_high = mean(margin_safety_high, 'omitnan') / std(margin_safety_high, 'omitnan'); % Calculate convergence, failure probabilities, and reliability indices for i = 1:length(years) for j = 1:length(scenarios) adjusted_water_level = total_water_level + slr_scenarios(i,j); % Failure probability failures = sum(adjusted_water_level >= critical_height_high); failure_prob_total_high(i,j) = failures / num_repetitions; % Reliability index margin_safety = critical_height_high - adjusted_water_level; mean_margin = mean(margin_safety, 'omitnan'); std_margin = std(margin_safety, 'omitnan'); reliability_indices_high(i,j) = mean_margin / std_margin; 102 % Convergence for k = 1:length(sample_sizes) n = sample_sizes(k); failures_n = sum(adjusted_water_level(1:n) >= critical_height_high); convergence_total_high(k,i,j) = failures_n / n; end % Find convergence point for k = window_size:length(sample_sizes) window = convergence_total_high((k-window_size+1):k,i,j); if std(window) < tolerance convergence_points_high(i,j) = sample_sizes(kwindow_size+1); break; end end end end % Save essential variables for height calculations save('base_monte_carlo.mat', ... 'random_hours', % For consistent time sampling 'astronomical_tide', % Base astronomical tide calculations 'residue_values', % Meteorological residuals 'Hs_initial', % Initial wave heights 'Hs_adjusted', % Wave heights after 3m threshold 'adjust_indices', % Indices for adjusted heights 'P_direction', % Random values for direction assignment 'months', % Month assignment for seasonality 'P', % Original random values for initial Hs 'P_adjusted'); % Random values for adjusted Hs 103 2. Matlab used functions 1. `clear all` - Removes all variables from the workspace memory. - Cleans up the environment before running new code. 2. `close all` - Closes all figure windows. - Used to ensure a clean visualization environment. 3. `load()` - Loads data from a MAT-file into the workspace. - Used to load 'my_trabucador.mat' in this script. 4. `randi()` - Generates random integers from a uniform distribution. - Used to generate random hour selections within total_hours range. 5. `cell2mat()` - Converts a cell array into a regular array. - Used to convert harmonic components data into numerical arrays. 6. `deg2rad()` - Converts angles from degrees to radians. - Used to convert phase angles for astronomical tide calculations. 7. `randn()` - Generates normally distributed random numbers. - Used in meteorological residue calculation. 8. `rand()` - Generates uniformly distributed random numbers. - Used for probability calculations in wave height analysis. 9. `deal()` - Assigns input values to output variables. - Used to assign seasonal constants (A, B, C). 10. `ismember()` - Returns array elements that are members of a set. - Used for checking valid wave directions and months. 11. `sum()` - Calculates the sum of array elements. - Used in various probability and counting calculations. 104 12. `cumsum()` - Calculates cumulative sum of array elements. - Used in wave direction probability calculations. 13. `discretize()` - Bins numeric data into categorical data. - Used in wave direction analysis. 14. `categorical()` - Converts data to categorical array type. - Used for wave direction categorization. 15. `sqrt()` - Calculates square root. - Used in various wave parameter calculations. 16. `zeros()` - Creates array of zeros. - Used to initialize various arrays throughout the script. 17. `nan()` - Creates array of NaN (Not a Number) values. - Used to initialize arrays for wave parameters. 18. `mean()` - Calculates average of array elements. - Used with 'omitnan' parameter for reliability calculations. 19. `std()` - Calculates standard deviation. - Used with 'omitnan' parameter for reliability calculations. 20. `fprintf()` - Outputs formatted text to display results. - Though not explicitly shown in the provided code snippet, it is commonly used for displaying results in MATLAB. 21. `floor()` - Rounds numbers down to the nearest integer. - Though not explicitly used in the provided code, it would be used for interval calculations. 22. `cell()` - Creates cell arrays. - Used in the script for handling the wave directions array - Example in code: `wave_directions = {'N', 'NE', 'E', 'SE', 'S', 'SW', 'W', 'NW'};`. 105 23. `length()` - Returns the length of arrays or vectors. - Used extensively in the code for loop iterations and array sizing. - Example in code: `for i = 1:length(years)`.