Coastal Flooding in the Balearic Islands During the Twenty-First Century Caused by Sea-Level Rise and Extreme Events

Sea-level rise induces a permanent loss of land with widespread ecological and economic impacts, most evident in urban and densely populated areas. Potential coastline retreat combined with waves and storm surges will result in more severe damages for coastal zones, especially over insular systems. In this paper, we quantify the effects of sea-level rise in terms of potential coastal flooding and potential beach erosion, along the coasts of the Balearic Islands (Western Mediterranean Sea), during the twenty-first century. We map projected flooded areas under two climate-change-driven mean sea-level rise scenarios (RCP4.5 and RCP8.5), together with the impact of an extreme event defined by the 100-year return level of joint storm surges and waves. We quantify shoreline retreat of sandy beaches forced by the sea-level rise (scenarios RCP4.5 and RCP8.5) and the continuous action of storm surges and waves (modeled by synthetic time series). We estimate touristic recreational services decrease of sandy beaches caused by the obtained shoreline retreat, in monetary terms. According to our calculations, permanent flooding by the end of our century will extend 7.8–27.7 km2 under the RCP4.5 scenario (mean sea-level rise between 32 and 80 cm by 2100), and up to 10.9–36.5 km2 under RCP8.5 (mean sea-level rise between 46 and 103 cm by 2100). Some beaches will lose more than 50% of their surface by the end of the century: 20–50% of them under RCP4.5 scenario and 25–60% under RCP8.5 one. Loss of touristic recreational services could represent a gross domestic product (GDP) loss up to 7.2% with respect to the 2019 GDP.


INTRODUCTION
Mean sea-level rise (MSLR) is one of the most certain consequences of human-driven climate change (Nicholls and Lowe, 2004). MSLR is quite relevant because of its potential impact over highly densely populated coastal zones, which also concentrate important natural and socioeconomic assets (Nicholls and Cazenave, 2010). It is expected that MSLR will partly submerge low-lying areas and increase coastal exposure to extreme events. According to projections, these changes will differ substantially among different regions due to the spatially varying mechanisms contributing to MSLR (such as ocean circulation and steric modifications, variations in the wind patterns, and water mass redistribution resulting from gravitational changes due to mass load variations) (Slangen et al., 2014).
The impacts of MSLR in the Mediterranean region are particularly worrisome since the population fraction living up to 10 m above mean sea level (MSL) reaches 34%, in contrast to 10% worldwide (Lionello et al., 2012). Moreover, the touristic boom experienced during the 1960s promoted an enormous population and urbanization boost on Mediterranean coastal areas, which continues nowadays. The phenomenon is more intense in sandy coastal environments, which provide the natural resource that attracts tourists to the region (Roig-Munar et al., 2019).
Beaches are among the most vulnerable ecosystems to MSLR (Vitousek et al., 2017;Vousdoukas et al., 2020), facing shoreline erosion and coastal flooding. Under natural conditions, beaches can adapt to MSLR by retreating, provided enough landward space and sediment supply are available (Cooper et al., 2020). Instead, most Mediterranean beaches lose width (i.e., backshore surface narrowing) due to the lack of accommodation space caused by heavy urbanization.
Beaches play an important role as natural coastal defenses, so their retreat and eventual disappearance increase the hazard vulnerability of the coastal region. Moreover, beach narrowing implies a loss in beach environmental services, which are critical to the economy of tourist destinations, since recreational services are dependent on the beach backshore functional surface for recreational activities (e.g., Valdemoro and Jiménez, 2006;Jiménez et al., 2007).
In this work, we assess the physical impacts of MSLR along the sandy coast of the Balearic Islands archipelago (Western Mediterranean sea), adopting a regional approach. Projections of MSLR, outputs of basin-scale wind waves and storm surge hindcasts, as well as local high-resolution topographic information, are combined to assess the long-term (up to 2100) flooding and erosion along 160 km of sandy beaches. Coastal flooding results are classified and quantified in terms of land use to aid the development of future adaptation mechanisms. The Balearic Islands are a well-known tourist destination dependent on sun-and-beach recreation activities. Furthermore, an economic assessment of the loss in recreational services is provided, based on shoreline retreat results. This assessment translates into monetary terms the effects of climate change over those beaches, which may lead to future exploration on how the Balearic Islands economy and society will be affected by climate change.
The paper is structured as follows. Section 2 presents the geographical context of the Balearic Islands, their coastal areas, and beaches. Section 3 describes the methodology used to obtain the necessary data and perform the flooding, erosion, and monetary analyzes. Sections 4, 5 present the obtained results. Finally, section 6 reviews the effects of coastal flooding over the different land uses, and section 7 quantifies the economic impact of beach erosion.

STUDY SITE
The Balearic Islands (Western Mediterranean sea) form an archipelago composed of four major islands: Mallorca, Menorca, Eivissa, and Formentera; and some small unpopulated islets. They constitute a prolongation of the Baetic mountain range, more than 100 km away from the east coast of the Iberian Peninsula and oriented in a SW-NE distribution. Eivissa and Formentera constitute a differentiated subarchipelago known as the Pitiüses (Figure 1).
The archipelago has about 1,700 km of coastline, where 867 beaches make up 10% of that length. The typical Balearic beach appears sparsely along the coast, accommodated in the space allowed by a physiographic control. It is short and narrow; presents features such as rocks, cliffs, and islets; and tends to be at the bottom of a wall-sided embayment. However, 2.4% of the Balearic beaches are longer than 1 km (Gómez-Pujol et al., 2019). According to Gómez-Pujol et al. (2019), who assessed regional shoreline change at the Balearic Islands from 2002 to 2012, 80% of beaches are stable, meaning they do not erode nor accrete more than 50 cm per year.
According to Wright and Short 1984 classification, approximately 20% of the Balearic Islands' beaches are enclosed and present reflective conditions; 42% of them are semi-enclosed, exhibiting intermediate but reflective-skewed conditions; 27% are exposed, non-protected beaches . The remaining 11% corresponds to anthropic or artificial beaches (generated by the presence of a groin or a dike).
Regarding beach sediment, most beaches are composed of sand, although more than 20% of them are mixed or biggersediment beaches. There is no substantial sediment fraction of fluvial origin since water streams in the Balearic Islands are ephemeral. Moreover, the connection to the mainland undergoes depths greater than 800 m, so sediment has a local origin. About 50-75% of sediment is bioclastic, although up to 10-30% of it can be of lithoclastic origin depending on the availability of cliff-detached material . The leading producer of biogenic material is the endemic seagrass Posidonia oceanica (Fornos and Ahr, 1997). It colonizes nearshore zones all around the Balearic Islands and provides housing to a variety of marine species. It also stabilizes submerged beach sediment and attenuates wave energy (Infantes et al., 2009).
Concerning marine forcings, the Balearic Islands' beaches are microtidal: the spring tidal range is smaller than 25 cm (Orfila et al., 2005). Regional winds are moderate and have a short fetch. As a result, they produce moderately short period waves: the typical significant wave height (H s ) range is 0.1-1 m and the associated periods are between 3 and 6 s (Alvarez-Ellacuria et al., 2011). The mean incoming wave direction is from the north and the north-west, caused by the Pyrenees and the Ebro valley (see Figure 1). The most energetic conditions occur between December and January, and the least ones during June and August (Cañellas et al., 2007). However, frontogenetic activity is very relevant in the basin (i.e., many atmospheric fronts are generated in the region), producing a large variability in wave regime (Morales-Márquez et al., 2020). Beach morphodynamics is mainly controlled by wave climate: the cross-shore coordinate is controlled by wave height and wave period, while the alongshore also depends on wave energy flux. For the few cases with available data, beach planform is stable (meaning no significant changes occur in the wave energy flux direction) and presents a seasonal cross-shore variation, which consists of an aerial beach loss during severe storms and a gradual recovery during mild wave conditions (Enríquez et al., 2017;Morales-Márquez et al., 2018).
Beach tourism has a large impact on the economy of the Balearic Islands, accounting for a 35% of the total Gross Domestic Product (GDP) (Llompart et al., 2020). In 2019 the Balearic Islands received 16.4 million visitors, mostly seeking beach and recreation during summer. This quantity represents more than 10 times the amount of archipelago inhabitants and makes the Balearic Islands one of the most popular destinations in Europe (Llompart et al., 2020).
Regional coastal management has focused on providing leisure services and keeping the backshore functional for recreational uses, with activities such as terrain smoothing, mechanical cleaning, Posidonia oceanica beach wrack (Balestri et al., 2006) withdrawal, and sand nourishments (Roig-Munar et al., 2019).
In the following, we analyze beach erosion and coastal flooding associated with climate change over 192 beaches

Regional Mean Sea Level, Storm Surges and Wind-Waves
The marine forcing around the Balearic Islands has been characterized under different scenarios at the regional scale. For flooding analysis, the magnitude of interest is the total water level, defined here as: where MSL is the mean sea level projected for the different scenarios considered, WS the component associated with the wave setup, and SS the storm surge. The latter two consist of 3-h time series, which add higher-frequency variability. Astronomic tides were omitted in this study. The portion of wave run-up not associated with wave setup was not considered because we focus FIGURE 2 | MSLR projections under RCP4.5 (A) and RCP8.5 (B) climate scenarios, calculated according to Kopp et al. (2014). Black lines indicate the multi-model ensemble median and shadowed regions indicate the 17-83% and the 5-95% probability intervals. The four colored lines indicate the MSLR evolutions considered for the analysis of beach erosion, while the six values labeled in bold indicate the MSLR cases considered for the analysis of coastal flooding.
on coastal extreme water levels that cause temporary flooding at time scales of hours to days. On the contrary, wave run-up acts on shorter time scales around a mean value of sea level. MSL projections were obtained following Kopp et al. (2014). Figure 2 represents local MSLR projections for RCP4.5 and RCP8.5 climate scenarios. We considered the multi-model ensemble 17 and 83% probabilities, defining a central 67% probability interval ("likely" range in the IPCC report of Stocker et al., 2014).
For WS and SS, we used data from Mentaschi et al. (2017) to Vousdoukas et al. (2017), who computed waves with the windwave spectral model WaveWatch III (Tolman, 2002), and storm surges with the hydrodynamic model Delft3d-FLOW (Deltares, 2006), consistently forced by atmospheric pressure and surface wind fields from ERA-Interim reanalysis (hindcast spanning the period 1979-2014) and from 6 CMIP5 GCMs for the historical period  and future projections (2070-2099). The temporal sampling of the hindcast is 3 h for H s and 6 h for SS. Figure 3 shows the dynamic models' grid points around the Balearic Islands, depicting the differences in the 50-year return period of H s and SS between the projections and the historical records. Multi-model mean differences for H s are shown in panels Figure 3A  Results indicate that the dispersion of extreme values among the models is larger than the expected changes in most locations. Moreover, it is unlikely to extract a realistic uncertainty from such a limited number of models (the datasets we use for H s and SS, from Vousdoukas et al. (2017), are obtained from 6 climate models: ACCESS1.0, ACCESS1.3, CSIRO-Mk3.6.0, EC-EARTH, GFDL-ESM2G, GFDL-ESM2M).
Consequently, given the unclear tendency and magnitude of the potential change plus the lack of precision in the uncertainty characterization, we decided not to account for any projected change in sea-level extremes. In other words, we assume that wave climate and storm surges will remain unaltered in the future. Therefore, we only used the hindcast outputs for both waves and storm surges, even to compute extremes. This way, we focus on the largest uncertainty (that from mean sea-level projections) and do not consider the others, such as those from storm surges and waves, which is in agreement with a fraction of the existing studies (e.g., Toimil et al., 2017b;Sanuy et al., 2018).

Wave Propagation
We selected 11 points (virtual buoys) that face all islands' orientations from the wave hindcast described above (see Figure 4). Deepwater conditions of those virtual buoys were propagated up to a set of reference points located about 30 m in depth, considered as representative of nearshore conditions. A map with the reference points locations (indicated as purple stars) is provided as Supplementary Material.
Due to the high computational cost of the 3-h 36-year-long wave hindcast numerical propagation, linear wave theory was applied to build time series of nearshore waves, as described in detail in the following.
Considering a wave propagating over a slowly varying bathymetry, before breaking and in the absence of wind and diffraction, its wave height H r at depth h r (before breaking) is related to its deepwater wave height H 0 as Dean and Dalrymple (1991): where K s is the shoaling coefficient (which depends on local depth and wave period), and K r is the refraction coefficient Frontiers in Marine Science | www.frontiersin.org (which depends on the bathymetry and the direction of the wave at deep waters, θ 0 ). The shoaling coefficient relating the variation in wave height from deep to intermediate waters is given by: where g is the acceleration of gravity, T is the wave period (which is not modified during propagation), h r is the water depth at the propagation destination point, k r is the wave wavelength, and ξ r = 2k r h r / sinh 2k r h r . In order to solve Equation (2), it is necessary to infer K r (θ 0 ) and, since this equation relates the change in wave height between two specific (though arbitrary) locations, we need to do so for each pair of virtual buoy and reference point. The process to obtain K r is explained in the following paragraphs. Hereafter, we consider all wave heights as significant wave heights (H s ), all periods as peak periods, and all directions as peak directions. We classified all hindcast sea states of each deepwater virtual buoy (red dots in Figure 4) in eight octants (waves coming from the N, NE, E, SE, S, SW, W, and NW). Each octant's monthly maximum wave height was propagated (with their corresponding period and direction) over 500 m resolution meshes (shown in Equation 4) using the action balance model SWAN (Booij et al., 1999), by forcing the described sea states along all the mesh side of that virtual buoy and neglecting the effects of wind. Propagations were conducted independently for Mallorca, Menorca, and the Pitiüses, defining a total of 42,528 propagations (3 domains × 4 deepwater points × 8 octants × 443 months).
An example of the SWAN configuration file, the code used to generate SWAN forcing files, and the code used to order all SWAN runs are included as Supplementary Material.
For each reference point and deep water virtual buoy, we computed the ratio K r = H r /H 0 · 1/K s relating the wave height at the virtual buoy and the reference point. K r (θ 0 ) was estimated by polynomial fitting, using the computed ratios for wave height of all octants at the same time. After a trial and error process, we found a fitting of the form: with p 1 , p 2 , p 3 , and p 4 parameters to be adjusted, abs the absolute value and θ 0 (the virtual buoy peak direction) between -180 • and +180 • . Figure 5 shows an example of fitting, which is representative of the majority of wave propagations. Four time series of wave height for each reference point were obtained applying Equation (2) and the fitted K r (θ 0 ). The final 3-h H s time series, representative of beach wave conditions, was built as the timewise maximum of its four reference points H s time series. The corresponding T series was defined as the peak period associated with the selected H s at each time. In line with the regional approach, wave trains were assumed to reach the coast perpendicularly. For this reason, nearshore wave direction time series were not computed. Also, to facilitate the understanding of this explanation, a graphical summary of the wave propagation process is presented as Supplementary Material. Many studied beaches receive wave trains from multiple deepwater points (especially those oriented toward the computational FIGURE 4 | Computational meshes used for wave propagation. Red dots indicate wave hindcast data points . Bigger dots indicate the data points we used for propagation (virtual buoys). Meshes are colored according to the bathymetry. domain vertices). However, wave trains that impact an island on one side should not affect the opposite side. We could manually decide which deep-water points should be considered for each beach, but this procedure is subjective and arbitrary. Instead, we assumed that wave trains with any direction coming from any deep-water point could reach (potentially) all beaches and thus we computed all the associated K r fittings. For the cases where the bulk of and island blocks wave train propagation, the obtained K r fitting presents small values (near to zero), and the resulting nearshore wave heights associated with those wave trains (computed with Equation 2) are negligible. Since the final nearshore wave height is the maximum among the four propagated waves from deep-water points, our model accounts for shadowing intrinsically.
Regarding diffraction, SWAN does not solve diffraction directly but uses a phase-decoupled approach. This method provides the qualitative behavior of diffraction processes for a restricted set of conditions listed in the technical manual of SWAN (the SWAN team, 2019a). In principle, features large enough to be represented in the 500 m computational meshes induce diffraction in the SWAN model (which can be captured by Equation (2), since K r is adjusted considering SWAN outputs). However, SWAN needs a grid size five to ten times smaller than the wavelength of a wave train to compute its diffraction adequately the SWAN team, 2019b). On the contrary, smaller features are not represented in the computational meshes, and thus their diffraction can not be estimated. Unfortunately, high resolution bathymetries are only available for specific locations. Considering all these reasons, we did not use the diffraction option of SWAN (DIFFRAC). Handling of diffractive obstacles is out of the scope of our regional approximation and is left to future, local studies.
In order to validate the reference points' H s , we compared our time series with the available in situ data at three coastal locations. The root mean square error at these locations ranges 0.20-0.28 m for H s and 1.04-1.54 s for T, while the bias (generated time series with respect to measures) takes values between -0.16 and -0.05 m for H s and between -0.29 and 0.01 s for T. The results and details of the comparison process, as well as the definition of the metrics used, are contained in Appendix A.
FIGURE 5 | Black dots represent refraction coefficients K r = H r /H 0 · 1/K s computed for one pair of deep water point (virtual buoy) and one reference point vs. virtual buoy peak wave direction (θ 0 ) for each case. Magenta line indicates the polynomial fitting [Equation (4)] performed over those points.

Computation of the Coastal Flooding Extent
Flooding is conducted over multiple zones, selected by considering the beaches under study and the topography, using terrain with high elevation as boundaries. A map with the location of all the flooding zones is provided as Supplementary Material. We employed the topography of "Servei d'Informació Territorial de les Illes Balears" (SITIBSA), as given by its 2 m resolution Digital Surface Model (DSM). The DSM is derived from LiDAR data of 0.5 points/m 2 point density and whose vertical RMSE is at most 20 cm. It can be consulted at (https://www.caib.es/sites/sitibsa/es/n/mdt-70457/).

Permanent Flooding Induced by MSLR
The values of MSLR used to quantify the permanently flooded coastal areas in 2050 and 2100 were the 17th and 83rd percentiles of both RCP4.5 and RCP8.5 multi-model ensemble spreads (Figure 2).
Permanent flooding was modeled using a bathtub method. Classical bathtub models consider all areas with an elevation below a given MSL as flooded (Yunus et al., 2016;Teng et al., 2017), being computationally cheap and thus suitable for regional assessments. However, they usually neglect physical processes, hence reporting greater flooding extents (Seenath et al., 2016). Aiming to diminish those undesired effects, we performed a multi-step correction over the bathtub results.
The application of our bathtub approach is as follows. First, we use the DSM to create a binary mask under each MSL value, where each DSM pixel is classified as flooded or not flooded depending on whether its elevation is lower or higher than the MSL considered, respectively. This binary map is then manually corrected to ensure surface hydrological connectivity, characterizing those elements misrepresented in the DSM (irrigation channels, thin walls, etc.). Afterward, we identify which binary flooding mask pixels are connected to the sea through a flood-filling algorithm. Flood-filling algorithms take an input point and find all pixels presenting the same value of a certain variable (what is called the region of that pixel), surrounded exclusively by pixels that present a different value on that variable (Rogers, 1998). Here, the initial point is an arbitrary pixel over the open sea, while the region comprises those flooded pixels connected to the initial one by a path made entirely of flooded pixels. After flood-filling, all pixels outside the region are marked as not flooded. Finally, small isolated groups of nonflooded pixels were also classified as flooded since they usually indicate undesired effects over the elevation map. Specifically, we reclassified groups of ≤50 pixels (≤200 m 2 ) if their elevation was lower than 20 cm above the MSL considered, and groups o of ≤100 pixels (≤400 m 2 ) otherwise. For each MSL case, we defined a new coastline as all non-flooded pixels connected with flooded pixels. This coastline was the one considered during the temporary flooding simulations. The whole modified bathtub process is summarized in the Supplementary Material.

Temporary Flooding Induced by Coastal Extreme Sea Levels and MSLR
Coastal extreme sea levels arise from the combination of storm surges (SS) and wave setup (WS), i.e., WS + SS, over the corresponding MSL. The spatial extent of coastal flooding induced by coastal extreme sea levels was simulated, at each flooding zone, with the LISFLOOD-FP model (version 7.0.6, http://www.bristol.ac.uk/geography/research/hydrology/ models/lisflood/), which is based on the shallow-water equations (Vreugdenhil, 1994). We performed 1-day simulations over the DSM described above, using the acceleration mode (i.e., neglecting convective acceleration but keeping the local acceleration term in the shallow-water equations), a constant Manning parameter of 0.06, and the default diffusion controlling weighting factor of 0.9 (which modifies the discharge time derivative to incorporate information about the discharge in surrounding cells in such a way that diffusive effects appear, thus reducing numerical instabilities) (de Almeida et al., 2012). An instantaneous sea-level time series was used to force the model along the coastline defined for every MSL value (the ones resulting from our permanent flooding estimation). To apply LISFLOOD-FP over the totality of each flooding zone (this is, including the regions between the beaches of the flooding zone), a single WS + SS time series was forced within each flooding zone, defined as the weighted mean of WS + SS from all the beaches inside this flooding zone, where the weights are the alongshore sizes of each beach. Following Stockdon et al. (2006), the WS of each beach can be estimated as: with β f being the foreshore slope of the beach, and L 0 the peak wavelength at deep waters. The value 0.35 is a constant calibrated using experimental data. We remark that for the second equality of Equation (5) we have assumed deep water conditions. Here, we propose a modification of this equation as: i.e., incorporating reference point's significant wave height (H r ), where α is a factor to be calibrated. H r is used instead of H 0 to remove those deepwater conditions that do not reach the coast due to their propagation direction. Factor α is obtained by dividing Equation (5) by Equation (6) and using a characteristic value for the ratio H r /H 0 according to the results of the numerical wave propagation (the value corresponding to the average angle of the H r /H 0 slopes bigger than 0.05). Foreshore slope β f is estimated using the distance between the coastline and the -10 m isobath, obtained from the nautical charts of the "Instituto Hidrográfico de la Marina" (this isobath is the shallowest one reported for the set of studied beaches in those nautical charts). We used a different β f value for each beach, and within each beach we used the same value for all the studied cases (i.e., we assumed that β f does not change with time or with MSLR).
Extreme values from the weighted WS + SS time series were characterized by fitting a Generalized Pareto Distribution to its exceedances over the 99th percentile. Independence among events was ensured by declustering with a time window of 72 h. The magnitude of the extreme water level event was defined as the fitted Generalized Pareto Distribution 100-year return level and its duration as the median duration of the exceedances (events exceeding the 99th percentile).
Finally, the time series forced over the flooding zone coastline was constructed as a triangular pulse over the corresponding MSL with an intensity and a duration as defined above (we keep sea level at the corresponding MSL after the end of the triangular pulse, until the simulation ends). The decision to use a triangular shape for the forcing is justified in Appendix B. An example of the LISFLOOD configuration for one of the flooding zones is provided as Supplementary Material.
Note that our LISFLOOD setup was not validated in the flooding zones studied due to the lack of measurements. Instead, we performed a sensitivity analysis in one flooding zone (the one presented in of the Supplementary Material). In our case, the only free parameter is the Manning coefficient (since all the other inputs are study design decisions that have been already discussed: mean sea level; extreme event shape, duration, and intensity). The Manning coefficient was varied from 0.006 to 0.6 with a change in the flooding extent lower than 0.04% for MSLR=0; 0.1% for MSLR=18 cm; 5.9% for MSLR=32 cm; 6.9% for MSLR=36 cm; 23.6% for MSLR=46 cm; 4.6% for MSLR=80 cm; and 0.36% for MSLR=103 cm (all of them with respect to the original simulations with Manning coefficient of 0.06). Changes are small except for the case of MSLR=46 cm, coinciding with a big area presenting an elevation near to this value and thus being more sensitive to the Manning coefficient.
Total flooded areas reported in this manuscript are the sum of the flooding extents categorized as one of the four land-uses we considered. Inner water bodies are not well characterized in the DSM because they are not identified as such. Moreover, interpolation errors occur over them. By removing the areas marked as LU-WB in the Corine database those effects are reduced, since those areas currently flooded (permanently or periodically) are not included in the results.

Computation of Beach Erosion
Balearic Islands beach sediment has a biogenic origin, without river income. Thus, all the sediment involved in beach erosion and accretion is already within the system. Under these conditions, morphodynamic variability can be split in the cross-shore and alongshore directions. Beaches accommodate maximizing wave energy dissipation, rotating until its alongshore direction is normal to the mean wave energy flux at the breaking point. Cross-shore variability (i. e., beach profiles) is the response to waves (shorter time scales) and to MSL (longer time scales).
That being said, the most accurate way to evaluate a specific beach is to perform a whole 3D simulation taking into account all changes inherent to the different scenarios. However, that is not affordable, neither for a specific beach (up to the 2100 horizon) nor for the regional approach we adopted. Therefore, we were forced to assume some simplifications to the problem.
In the first place, as stated in section 2, only 27% of the Balearic Islands' beaches are open (and thus can freely exchange sediment with its surroundings). Accordingly, we neglected the alongshore transport between the beach and its surroundings. In the second place, generalized information about wave energy flux does not exist, so beach planform changes can not be assessed; moreover, beach planform is stable for the few cases where this information is available. Therefore, alongshore averaged cross-shore beach width was used. In the third place, the number of recreational activities occurring at beaches depends directly on the aerial beach surface, equivalent through the alongshore beach size to the average cross-shore beach width. Since we focused on an alongshore-averaged magnitude, internal alongshore transport (i.e., occurring between the beach and itself).
Because of these reasons, we quantified beach erosion as the average cross-shore beach width recession. Classically, the Bruun Rule (Bruun, 1962) has been used for this purpose. It indicates that the beach tends to an equilibrium state whose changes in shoreline can be described (in the form presented by Dean, 1991) as: where y is the (seawards) change in equilibrium beach shoreline position with respect to a reference sea level, S indicates the elevation of sea level with respect to that reference, B is the berm height, h * is the breaking depth, and W * is the active profile width, computed as W * = (h * /A) 3/2 (Dean, 1991), where A = 0.067 ω 0.44 s , with ω s being the settling velocity of mean sand diameter (D50) (van Rijn, 1984). According to similarity theory, the following relation holds at the breaking point: H b /h * = κ, where H b the wave height over breaking depth, and κ is the breaker index (which we assumed to be 0.71).
As stated in section 2, most of the Balearic Islands' beaches lack accommodation space, either because of heavy urbanization or because of natural physiographic controls (pocket beaches in front of cliff walls). In line with our regional approximation, we assumed the backshore of all beaches would remain fixed, so shoreline position changes equal to cross-shore beach width changes ( y = w). Taking this into consideration, using the current mean sea level (MSL 0 ) as the sea-level reference, and neglecting changes in sea level other than MSLR, Equation (7) can be rewritten as: w 0 is the equilibrium beach width for the case MSL = MSL 0 and w the equilibium beach width for a MSLR equal to MSL − MSL 0 . Bruun's Rule has some limitations since it assumes a constant wave climate, a constant sediment budget, and the existence of a depth of closure, and also neglects the effect of longshore sediment transport gradients. However, given the nature, dynamics and physiography of the analyzed beaches, as well as the regional context of the study, Bruun's Rule provides an acceptable first approach to the expected beach width changes.
Regarding to the rest of variables needed, according to Gómez-Pujol et al. (2019) beaches in the Balearic Islands are in equilibrium, and so we assumed w 0 as the one measured by the most recent orthographic aerial photographs. For each beach, we estimated a single B value considering the LiDAR point cloud used to generate SITIBSA's DSM. Wave height was computed at different depths by propagating the wave height time series of the beach reference point (H r ), first to deep waters and then to the desired depth. We accomplished that using Equation (2) and assuming K r equal to one.
Following Toimil et al. (2017a), we applied Equation (8) to a set of multiple synthetic time series of H b , thus producing multiple potential realizations of beach width evolution. Our synthetic H b time series are statistically consistent with the corresponding reference point time series at every site (i. e., they have the same statistics but with a different chronology). Since the set of synthetic time series generated is representative of the potential future conditions, the resulting realizations of beach width time series can be used to estimate the statistics of the future beach width evolution. The details of the synthetic time series generation process and its validation are described in Appendix C.
For each beach, we generated 500 synthetic H b time series and inputted them into Equation (8), in combination with the projected evolutions of MSLR corresponding to the 17 and 83% of RCP4.5 and RCP8.5, obtaining 500 time series of w for each MSLR evolution, from which mean beach width evolution and standard deviation were computed. Since both metrics presented a strong seasonality, resulting time series were filtered with a 5-year running average.

Beach Recreational Services Monetary Valuation
Beach narrowing due to erosion resulting from our analyzes was used to estimate the decrease or loss of beach recreational services in monetary terms. We adopted the results presented in Enríquez and Bujosa Bestard (2020), who measured the economic impact of beach loss on beach tourism through choice experiments focused on the Balearic Islands. Other works also assessed the recreational services provided by beaches at the Balearic Islands (Pérez-López and Roig-Munar, 2007;Riera et al., 2007), but their surveys were not focused on climate change impacts or on the tourists that comprise the largest portion of yearly beach users in the entire archipelago. Enríquez and Bujosa Bestard (2020) specifically assessed tourists' willingness to pay for the introduction of policies aimed at reducing climate change impacts.
Choice experiment methodology is a stated preference nonmarket valuation methodology based on experimental surveys. In these surveys respondents are presented with different scenarios describing specific changes in the levels of the good under valuation (in this case the surface of aerial beach) and are asked to choose the scenario providing the highest level of wellbeing. This is a standard method to assess economic impacts of climate change (Shoyama et al., 2013;Andreopoulos et al., 2015;Remoundou et al., 2015;Rodrigues et al., 2016;Chaikaew et al., 2017;Torres et al., 2017;Faccioli et al., 2019).
In their assessment, they estimated that tourists are willing to pay 1.23 e/day of stay (on average) to recover one meter of beach shoreline to compensate beach retreat due to MSLR. Using this value we estimated the loss of beach recreational value in monetary terms. In this regard, we combined the data obtained in section 3.4 with the Enríquez and Bujosa Bestard (2020) results, as well as with different attributes of each beach (beach location, type, accessibility, etc.) to determine the beach recreational economic value loss (BREV) as follows: BREV n e yr = WEBS n m 2 n WEBS n m 2 · TOUR n (tourists) ·WTPV e m · yr · ER n (m) where BREV n is the loss of beach recreational economic value for a specific beach, WEBS n is the weighted effective beach surface for a specific beach, TOUR n is the number of tourists (corresponding to the 2019 official estimate of the Balearic Islands Institute of Statistics, IBESTAT), WTPV is the amount of money each tourist is willing to pay in constant monetary units for an individual day of stay considering the average stay length of 6 days (Llompart et al., 2020), and ER n is the averaged shoreline retreat for each beach obtained in section 3.4.
Beaches in the Balearic Islands tend to be high-frequented (Mas Parera and Blázquez-Salom, 2005;Roig-Munar et al., 2020). In a small island (as each one of the Balearics, with travel times from one side to another lower than 1 h) tourists move from one beach to another. For this reason, the number of users is mainly related to the available dry beach surface. The articles cited above also indicate that the most frequented beaches are those that are in front of tourist stations and also those iconic virgin beaches publicized in touristic guides and brochures. Therefore, according to this descriptive data, we obtain the weighted effective beach surface for a specific beach by means of: WEBS n m 2 = S n m 2 · F type · F access (10) where S n is the beach surface for a specific beach; F type is a correction factor that introduces the typology of beach that, ranging from 0.5 for isolated beaches to 1 for urban beaches and isolated iconic beaches (Table 1), as the two latter ones comprise most of the demand (Roig-Munar et al., 2020); F access is a correction factor related with the beach accessibility (access, public transport, parking facilities) and beach recreation services (beach cleaning, sun huts, lifeguards...), creating a variation from 0.05 to 1 according to Table 1.

RESULTS FOR COASTAL FLOODING
The results of areas permanently and temporarily flooded due to MSLR combined with storm surges and extreme waves are summarized in Table 2, separately for each region, and classified by its land use. Note that temporary flooding due to extreme events is indicated as increases with respect to permanent flooding induced by MSLR only. Total extent of flooding in the case of extreme events is the addition of the two. Also, note that results are listed in terms of increases in MSL rather than in time horizons. All values reported in the table correspond to either mid or late twenty-first century (see Figure 2 for the equivalence between timing and value of MSLR). The complete set of flooding extent masks can be consulted online at https:// ideib.caib.es/impactes_costa_canvi_climatic/. Finally, note that those are estimates for the case that no action is taken. Figure 6 shows the permanent flooding extent for a MSLR of 103 cm over the entire archipelago, classified by municipalities. The area permanently flooded increases linearly with MSLR. Namely, in Mallorca flooding extent changes 24 ha/cm of MSLR (R 2 of 0.977), in Menorca 2 ha/cm of MSLR (R 2 of 0.975) and at  the Pitiüses the rate reaches 33 ha/cm of MSLR (R 2 of 0.956). Due to the different coastal lengths, absolute values of flooded extents are clearly dissimilar, ranging between areas up to thousands of hectares in Mallorca to a few hundreds of hectares in Menorca and Pitiüses, for MSLR above 0.5 m. Temporary flooding caused by coastal extreme sea levels increases mostly in Mallorca and Pitiüses (100-300 ha) and is smaller in Menorca, due to the more elevated coastal land. In terms of relative impacts of flooding, by mid-century and under RCP8.5 (MSL between 18 and 36 cm), permanent flooding reaches between 4.9 and 8.6 km 2 (0.10-0.17% of the archipelago surface), while extreme flooding adds an extra loss of 0.06 -0.07% to the archipelago land surface. This represents an increase of 16-53% with respect to the flooding extent of an extreme event occurring today. We expect these quantities to be very similar to those of RCP4.5 for the same time horizon.
Likewise, by the end of the century, under the RCP4.5 scenario (MSLR between 32 cm and 80 cm), the permanent flooding extent would fall between 7.8 and 27.7 km 2 (0.16-0.56% of the archipelago surface). Extreme flooding represents an extra 0.07-0.10% of archipelago surface, which is 44-96% bigger than the extreme flooding extent of the present. Under RCP8.5 scenario (46 cm and 103 cm limits) permanent flooding can extend up to 10.9-36.5 km 2 (0.20-0.70% of the archipelago surface) and extreme flooding adds 0.11-0.12%, an increase of 105-141% with respect to the extreme flooding that may occur nowadays.

RESULTS FOR BEACH EROSION
A compilation of the results for all the beaches is represented in Figure 7 as the CDF of the relative widths with respect to their present state for mid and late twenty-first century. Note that those are estimates for the case that no action is taken.
Around 2050, beach width distributions obtained for the lower and upper limits of RCP4.5 and RCP8.5 are similar. In Mallorca and Menorca around 3-15% of beaches are projected to lose half of their current surface, 0% in the Pitiüses. Between 25 and 55% of the beaches are projected to reduce their present area by less than 10%. Changes are significantly larger by the end of the century. Under the RCP4.5 scenario up to 10% of the Balearic Islands' beaches will lose 90% of the current area. 20-50% of the beaches in Mallorca, 15-50% in Menorca and 5-40% in the Pitiüses are projected to reduce their current size by half. In the worst-case scenario, under RCP8.5 by 2100, up to 15% of the archipelago beaches are expected to lose more than 90% of their surface, while 25-60% of beaches in Mallorca, 35-65% of beaches in Menorca and 23 -60% of beaches in the Pitiüses will lose more than 50% of their current surface. In order to illustrate the spatial distribution of beach erosion, Figure 8 shows the beach surface loss between the present and the 2100 RCP8.5 upper limit of MSLR projection classified by municipality.

COASTAL FLOODING IMPACTS ON LAND USE
The classification of the extent of flooded areas by their land use, according to Corine 2018, is listed in Table 2. Note that those are estimates for the case that no action is taken. Our analyzes demonstrate that wetlands are disproportionally affected by coastal flooding in the Balearic Islands, being clearly dominant in Mallorca and Pitiüses. Overall, these areas represent 65-85% of the total permanently flooded surface. In Menorca, semi-natural areas and forests are the most impacted. Flooding of coastal wetlands increases fast as MSL rises up to 0.5 m, due to their typical low elevation gradient. Above this value these impacts stabilize, indicating that the totality of coastal wetlands is lost. It must be mentioned that we have neglected the dynamic response of coastal wetlands, which may accrete in response to mean sea-level changes, depending on the rising rate and the space availability for sediment accumulation (Schuerch et al., 2018). The latter is rather limited in our case studies, though, but still our estimates should be considered as an upper bound for the given MSLR values.
The areas most affected by coastal flooding, coinciding with the low-lying regions, are formed essentially by lagoons and salt evaporation ponds. There are five flooding zones that account for 85% of the extent of total permanent flooding in the worstcase scenario (MSLR of 103 cm). These are concentrated in the North of Mallorca Island (Alcúdia, Pollença), the South of Mallorca (es Trenc), and parts of Eivissa and Formentera. Therefore, coastal flooding impacts are unevenly distributed around the archipelago. These results are relevant for public policies regarding the prioritization of most vulnerable coastal areas to MSLR.
Frequently, large areas are projected to be flooded because of the presence of a narrow connection to the sea, such as channels, torrents, or salt work entrances. This is illustrated in FIGURE 8 | Aerial beach surface loss by 2100 under the RCP8.5 upper limit of MSLR, relative to that of today, classified by municipality. Beaches considered are also shown. Figures 9A,B, where Port of Alcúdia flooding zone is the most paradigmatic example: the vast majority of its streets and multiple channels allow flooding of the urban nucleus and the wetland area due to MSLR. The city of Eivissa (Figures 9C,D) includes a good example of an artificially guided torrent.
In these maps, public facilities such as schools, sanitary buildings as well as energy infrastructure are marked to illustrate the utility of the assessment to identify critical infrastructures at risk. Major impacts on infrastructures and equipment affect streets, secondary roads, and tracking pathways. This is quite evident in Figures 9A,B, where urban streets at tourist stations and rural secondary roads communications become flooded, especially for MSLR larger than 40 cm. Municipalities are the administration in charge of the most affected communication and transport infrastructure. The power plant at the Alcúdia Bay ( Figure 9A) or the Eivissa airport are not affected by MSLR, neither by the associated extreme events, but they will be resting closer to the shoreline. Permanent and eventual floods do not affect any sanitary equipment (Figure 9), while some educational equipment may be damaged at Alcúdia bay (1 school and 1 institute, Figures 9A,B) and Eivissa (1 school, Figures 9C,D) if MSLR surpasses 46 cm. Figure 6 shows the fraction of permanent flooding extent associated to each land use considered, by municipalities, for the case of 103 cm of MSLR, for the entire archipelago.

LOSS OF RECREATIONAL VALUE DUE TO BEACH EROSION
The estimated losses of recreational value in the beaches of the Balearic Islands are shown in Table 3 for each island separately and for the entire archipelago. Note that those are estimates for the case that no action is taken. The loss ranges from 415 to 827 million euros under RCP4.5 for 2050, from 481 to 932 million euros under RCP8.5 for 2050, from 1,000 to 2,067 million euros under RCP4.5 for 2100, and from 1,421 to 2,428 million euros under RCP8.5 for 2100.
Under the RCP8.5 scenario, beach recreational services loss by 2050 could represent between 1.4 and 2.8% of the 2019 Balearic  Islands GDP, while by 2100 it represents between 4 and 7% of the 2019 Balearic Islands GDP. The Balearic Islands GDP is estimated at 33,800 million euros for 2019. It should be noticed that this economic value refers to services provided for tourism recreation, not including the local population. Otherwise, beaches provide other additional environmental services, such as coastal protection. In monetary terms, the most damaged beaches of the archipelago under the RCP8.5 scenario in 2100 correspond with the largest beaches located in the major coastal basins of the island: Alcúdia beach in the northern coast of Mallorca and "s'Arenal" in the southern coast of Mallorca, with 574 and 169 million euros of recreational services loss, respectively. In Eivissa, "Platja d'en Bossa" is the most affected beach, with a loss of recreational services equivalent to 147 million euros, whereas the rest of the beaches of the island remain one order of magnitude below this value. The impact on Menorca beaches, in terms of recreational economic value, is lower than in the rest of the islands. The largest impact belongs to Son Bou beach (southern Menorca) with a recreational service loss equivalent to 38 million euros for 2100 under the RCP8.5 scenario.

CONCLUDING REMARKS
In this work, we estimated the extent of coastal flooding and beach erosion along sandy coastlines in the Balearic Islands as a result of projected MSLR during the twenty-first century. These coastal impacts were quantified for two climate change scenarios, namely RCP4.5 and RCP8.5, combining regional projected MSLR and marine extreme events caused by waves and storm surges. Our analysis constitutes the first study of mean sea-level rise impacts in the Balearic Islands at the regional scale, comprising all sandy coasts. Due to the size of the study site, the analysis was conducted as a regional approach, meaning that we considered a series of methodological approximations. Most of these approximations are related to beach erosion. We assumed that shoreline retreat was equivalent to beach width loss, while those beaches without an urbanized backshore will just displace landwards, as stated by Cooper et al. (2020). Also, we used Bruun's Rule for the computation of beach shoreline evolution which, although provides an acceptable first approach to the expected beach width changes, possess many limitations. For this reason, we believe future studies using more advanced models, such as that of McCarroll et al. (2021), should be conducted.
The approach permitted to identify the potentially most affected areas in the region: the wetlands of Alcúdia, Pollença and es Trenc in Mallorca, and the saltworks of Eivissa (accounting for 85% of the 103 cm of MSLR permanent flooding); as well as the number of beaches undergoing a critical state of erosion: up to 60% of beaches in the region may erode more than a 50% of its current width for a MSLR of 103 cm. Thus, our results can be used to inform the regional and national administrators about the most critical zones in order to prioritize their actuation.
The estimated loss of recreational value for the considered beaches is estimated between 1.4 and 2.8% of 2019 Balearic Islands GDP by 2050 and between 4.2 and 7.2% by 2100, under the worst-case scenario. Importantly, these values do not include the value provided by other beach services such as its role in coastal protection.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
LG-P, MM, and AO conceived the idea of the study with the support of PL. PL, LG-P, MM, and AO developed the methodology. PL produced the results with the support of LG-P, MM, and AO. All authors analyzed the results and equally contributed to writing the manuscript. LG-P is also funded by the CGL2016-79246-P MINECO-AEI-UE-FEDER project. PL was supported by the FPU19/03081 grant from Ministerio de Ciencia, Innovación y Universidades of Spain.

A. Wave Propagation Validation
This appendix describes the instrumental data and the method used to validate the wave time series. The instrumental data consists of the available nearshore wave time series, which comprises data from three study sites: Cala Millor in the East of Mallorca, Platja de Palma in the Southwest of Mallorca, and Son Bou in the South of Menorca. Data was acquired using in situ Acoustic Wave and Current meters (AWACs) deployed at coordinates listed in Table A1 and 18 m depth, managed by SOCIB. The AWAC configurations can be consulted on the SOCIB webpage.
Pairs of simultaneous measured and synthetic H s and T as well as the quantiles (q n ) associated with the percentiles (0,1,2,...,99,100) of each time series are considered to compute the following metrics: where x symbolizes H s or T; N is the number of time series samples considered for the analysis; x m i , x m and σ x m refer to the i th sample, the mean and the standard deviation of the instrumental time series, respectively; while x s i , x s and σ x s refer to the i th sample, the mean and the standard deviation of the synthetic time series, respectively. Also, q s n and q m n are the quantiles of the percentile n, for synthetic and measured time series, respectively (both for H s and T); while q s and q m stand for the mean of the quantiles of the synthetic and measured time series, respectively. The results of the seven metrics described are listed in Table A1. Figure A1 shows scatter plots and QQ-plots for these data pairs.

B. Justification of the Triangular Flooding Forcing
In order to define the form for the forcing of the temporary flooding, the WS + SS exceedances of all the flooding zones were considered. A prior visual inspection indicated that the shape of the exceedances tends to be triangular-like, at least at the sampling rate used.

C. Synthetic Storm Surge and Significant Breaking Wave Height Time Series Generation
To characterize the statistics of H b time series, it is necessary to consider the intensity of their samples, as well as the time distribution and persistence of storms, which in this case are defined as the events where the intensity of the forcing exceeds some threshold during a certain time, called persistence of the event (Toimil et al., 2017a). Once the statistics of the forcings are known, it is possible to generate synthetic time series that are statistically consistent with the original one, but with a different chronology, in this case based on Soares et al. (1996). These can be used to generate an estimation of the beach width time evolutions probability distribution, by means of a Monte Carlo method (by using Equation (7). We characterize the logarithm of H b (clipped below -2) instead of H b itself. The forcing time series present annual stationarity. Thus, instead of characterizing the statistics of the whole time series, we classified their samples depending on their position inside the seasonal period and used the fragments of time series corresponding to each season (December, January, and February; March, April, and May; June, July, and August; September, October, and November).
For each group, we estimated its cumulative density function (CDF). For the logarithm of H b we used a piecewise function with two intervals, where the range below the 99.9 th quantile was computed by linear interpolation, and the interval above that threshold was computed as 1 − (1 − 0.999) · (1 − CDF GPD ), where CDF GPD is a Generalized Pareto distribution CDF, fitted by least squares.
Once the CDF of each group was estimated, we applied the inverse transformation method (Soares et al., 1996) to made each group's data follow a Gaussian distribution. Then, we computed and solved the Yule-Walker equations associated with each group's gaussianized data, obtaining an autoregressive filter for each group (Kay, 1993).
To obtain synthetic time series we generated Gaussian white noise and filtered it using the previously calibrated autoregressive filter (the output of the filter presents the same time correlation as the data used to calibrate the filter, though it follows a gaussian distribution). Afterward, the inverse transformation method was applied to obtain a time series presenting both the same time correlation and marginal distribution as the original data. Finally, in order to mimic the seasonality of the forcings, the filter and CDF of each season were switched periodically between the ones of the different groups. A comparison between the original and multiple synthesized time series for one beach is summarized in Figure A2 in the form of different statistics.