ROLE OF SPATIAL VARIABILITY OF SOIL RESISTANCE IN ALONGSHORE VARIABILITY OF COASTAL BARRIERS RESPONSE TO SUPERSTORM SURGES

Sand dunes and other natural coastal barriers (e.g. barrier islands) represent important components of the defense system against consequences of storm surges. However, in many coastal systems, major storm surges represent important drivers of coastal erosion. Increased extreme events potentially result in accelerated coastal erosion, coastal barrier breaching, and coastal flooding. The response of a barrier to a storm surge is often determined by mutual interaction among the driving hydrodynamics, the subsequent morphodynamics, and the local geology, including spatial variations of subaqueous bathymetry and subaerial topography. However, the effect of alongshore variability of soil properties on the alongshore varying response is not yet considered. Therefore, this study examines soil parameters that may affect coastal erosion during major storm surges. Moreover, it applies a novel extension of the numerical model XBeach that accounts for spatial variation of soil properties to an artificial dune system of spatially varying soil permeability. Results showed that variability of soil permeability alongshore the dune results in alongshore varying resistance to erosion so that breaches may occur at the locations of less resistance that are corresponding to locations of higher soil permeability. Outcomes of the numerical simulations proved also that reduced soil permeability represents a nature-based solution that increases the resilience of natural defense systems during major storm surges by mitigating rates of coastal erosion.


Introduction
Recent extreme storm surges induced for instance by (i) hurricanes Irma that hit Florida in 2017 (Xian et al. 2018) and Katrina affecting Louisiana in 2005(Fritz et al. 2007), (ii) the 2010 storm Xynthia that impacted south-western Europe (Bertin et al. 2012), and (iii) the 2013/2014 winter storms that impacted the UK coastline (Masselink et al. 2016) have proven that coastal areas are highly vulnerable environments to coastal floods. Globally, extreme storm surges tend to occur along 55% of the world's coastlines . Moreover, recent studies on climate change and its relation with coastal erosion (e.g. Vousdoukas et al. 2016;Walsh et al. 2016) have reported that climate change might lead to changes in frequency, intensity, spatial extent, duration, and timing of weather events, possibly resulting in unprecedented extreme events. As a result, coastlines and coastal barriers may experience severe impacts from coastal storms. Storm surges may substantially influence coastal dynamics, leading to accelerated coastal erosion, breaching of coastal barriers, and coastal flooding (Didier et al. 2019;Elsayed and Oumeraci 2016;Enríquez et al. 2019;Passeri et al. 2018;de Winter and Ruessink 2017).
The accurate determination of dimensions and locations of a possible barrier breach is a very important topic to predict accurate inland discharges and flood spreading for many applications such as coastal flood risk assessments and vertical intrusion of saltwater behind breached barriers (see e.g. Elsayed and Oumeraci, 2018 and references therein). On the other side, it is a real challenge that arises from the fact that the breaching process represents a mutual interaction between (i) the hydraulic loading (e.g. storm surges combined with high-frequency wind waves), (ii) the alongshore complexities of barriers topography and fronting bathymetry and (iii) the type and properties of the barrier's soil, which may also vary along the barrier. The former two factors are frequently addressed during studies relevant to response of beaches and coastal barriers to storm surges (e.g. Cohn et al. 2019;Elsayed 2017;van der Lugt et al. 2019;Mickey et al. 2018;Splinter et al. 2018), while the latter factor and its influence on the modeling of beach and dune erosion gained far less attention in past studies. In fact, recent studies (e.g. Foortse et al., 2019;Lemmens et al., 2016;Mohr et al., 2018Mohr et al., , 2021 have experimentally proven that a soil's tendency to erode depends also on soil properties, e.g. at first order its permeability. Nevertheless, substantially more research is still needed to quantify the order of magnitude of soil resistance to erosion, to understand the underlying processes within the soil's grain structure and to then better account for the spatial variability of soil resistance to erosion. This study is therefore aiming at (i) examining the effect of spatially varying soil properties on providing alongshore varying resistance to erosion (ii) extending the XBeach model (Roelvink et al. 2009) to account for the spatial variability of soil resistance to erosion, and finally (iii) applying the extended XBeach to a case study of spatially varying soil properties in order to numerically quantify effects of spatially varying soil resistance on coastal erosion and predictions of possible breaches under extreme storm surges. Next, a brief overview over the relevant literature is presented

Modeling spatially varying effect of soil properties on coastal erosion
Spatial variability of soil properties and their incorporation into modeling of coastal response and dune erosion was not addressed comprehensively so far in state-of-the-art studies. In this section, the modeling limitations that arose due to omitting mapping of spatial variability of soil properties in modeling of coastal erosion are highlighted. Besides, the recent efforts to account for soil properties and the extension of the XBeach model to account for spatially varying soil properties are presented.

Importance of accounting for soil properties in modeling of coastal erosion
Shoreline evolution and dune erosion during storm surges take place over temporal scales ranging from seconds (due to the individual waves) to hours (duo to the time-varying storm-tide levels) (Montaño et al. 2020). Flow velocities of individual waves, as they run up and hit coastal barrier elements during storm surges, are significant; they may reach approximately 3 m/s, particularly those of the turbulent uprush and gravity-induced backwash components of the swash motions (Masselink and Russell 2006;Pujara et al. 2020). Moreover, flow velocities during overwash and overflow conditions as well as through breaches are even larger, reaching approximately 10 m/s (Bisschop et al. 2010;Elsayed and Oumeraci 2016). During such very high flow velocities, sediment transport rates by prediction modeling tools, e.g. XBeach, are often overestimated because of employing conventional sediment pick-up functions based on the erosion of single grains (i.e. the grain-by-grain pick-up of Shields, 1936), which are only applicable for low flow velocities up to 1 m/s (Bisschop, 2018 and references therein). As a result, models based in such conventional sediment pick-up functions may fail in predicting values of sediment loads accurately (van Rijn et al. 2019). Thus, predictions of possible breaching sizes and dimensions using numerical models (e.g. XBeach) may still lack reliability (Elsayed and Oumeraci 2017;van Rijn et al. 2019;De Vet 2014).
In order to overcome the shortcomings of such conventional pick-up functions, McCall et al. (2010) introduced an artificial sediment transport limiter as shown in Figure 1 that limits the effective Shields parameter θ at 1.0 [-]. This limiter (known in XBeach by the keyword smax) serves as a proxy for unknown factors that control resistance to erosion during high flow velocity regimes by setting an upper boundary for the destabilizing depth-averaged sediment stirring velocity . This incident velocity for inception of sediment motion is formulated in XBeach in terms of the effective Shields parameter , median sediment size 50 , the relative density of sediment s ( / ), gravitational acceleration and non-dimensional friction coefficient as follows: This limiter means that is not allowed to exceed the range of 1.0 m/s in the calculations of sediment transport in XBeach. In this case, sediment transport has an upper bound that cannot be exceeded during higher flow velocities, irrespective of the physical magnitude of the flow field because the value of is restricted at the limiting value (i.e. ≤ 1). This artificial limiter has been applied in many erosion and breaching studies (e.g. Gharagozlou et al. 2020;Harley et al. 2016;Muller et al. 2016;Passeri et al. 2018b) often showing moderate prediction accuracies as compared to experimental or insitu observation data. Other studies (e.g. Elsayed and Oumeraci 2017;Terlouw 2013;De Vet 2014) reported that the artificial limiter approach may result in unexpected side effects and may, in some cases, lead to significantly worse predictions than without any limiter. This problem of overestimation provided a first motivation to this study as it revealed the essence of next steps in erosion modeling: to develop physically-based limiters of erosion rates based on studying the role of soil properties and the interlocking among sediment particles in hindering the erosion under high flow velocities.
As a tentative attempt to improve the overestimation problem using physically-based limiters, Elsayed and Oumeraci (2017) introduced the concept of soil resistance to erosion. It further attributed the overestimation of erosion to excess shear stresses required to overcome grain-stabilizations and initiate motion of interlocked grains (Agudo et al. 2014;Reid and Frostick 1984). Thus, soil resistance to erosion was represented by an amplification factor R that accounts for the additional critical shear stress required to initiate sediment motion during high flow velocity regimes ( ≥ 0). The soil resistance factor was integrated into the conventional sediment pick-up function of Soulsby-Van Rijn (van Rijn 2007a;b;Soulsby 1997 (2) where and are the bed and suspended load coefficients respectively (see e.g. Deltares (2015) for details), is the local water depth and is the critical stirring (threshold) flow velocity for sediment movement that is based on the Shields curve (Shields 1936). The critical depth-averaged velocity is calculated based on the critical Shields parameter , as in Eq 4: 2 = , . 50 ( −1) (4) As schematically shown in Figure 1, the introduced soil resistance parameter is a tuning parameter aimed at increasing the critical Shields parameter and hence the critical stirring velocity in Equations 2 and 3 by a value based on the overestimation magnitude. Though this tuning parameter reduces the mismatch between the effective and critical stirring velocities in Equations 2 and 3, it still does not provide an explanation of the forces amongst soil particles resulting in resistance and the role of soil properties in augmenting these bonding forces and the soil resistance to erosion. This is further explained by the dilatancy theory of Van Rhee (2010) as described in the following section.

Soil dilatancy and permeability effects on erosion
Van Rhee (2007,2010) introduced the dilatancy theory to explain a delay in the initiation of sediment motion during high flow velocities. He demonstrated that sediment movement under high flow velocities (i.e. >1.0 m/s) is no longer a grain-by-grain movement as described by Shields (1936); it is rather a sheet flow motion (hereafter dilatancy-reduced erosion regime) as shown in Figure 2. Therefore, a granular material moving in sheet flows needs to dilate (i.e. expand in volume) by a slight increase of the porosity from a value 0 to a value 1 as it is sheared. The volume change generally leads to a drop in pore pressure at the top of the sand bed (Bisschop et al. 2010). This pressure drop introduces a hydraulic gradient ( ) and hence an inflow of water that hinders the entrainment of sediment and provides resistance to erosion.
The modified critical Shields parameter after considering the excess in shear stresses required for soil dilation θ cr ,D can be calculated according to the model of Van Rhee (2010)  ) . θ cr,s = (1 + ). θ cr,s Van Rhee's expression (Equation 5) predicts that the threshold shear stress for fine sediments may be higher than the empirical Shields curve by a value that depends on the hydraulic gradient ( ); in turn, it further depends on the erosion velocity v e , the hydraulic permeability and porosity before and after dilatation (i.e. 0 and 1 ). The erosion velocity v e refers to temporal change in bed level ∆ /∆ which increases with higher flow velocities (Sethi 2015). The geometrical parameter A is a coefficient depending on sediment geometry (A=3/4 for single particles and 1.7 for a continuum; See Van Rhee (2010) for details). Van Rhee's equation provides therefore an evaluation of the soil resistance R [-] as a function of soil properties and erosion velocity. Soil resistance to erosion may thus be defined based on Equation 5 as the additional shear stress required to initiate sediment motion during the dilatancy erosion regime using the Shields curve as a reference line, as illustrated in Figure 1, where resistance vanishes (i.e. R = 0 when θ cr ,D = θ cr,s ). The resistance to erosion diminishes in case of high soil permeability or the flow velocity and hence the erosion velocity decreases (Elsayed et al. 2021). Therefore, the dilatancyreduced erosion regime starts to dominate at relatively high flow velocities (> 1.0 m/s), corresponding to higher bed shear stresses and higher ratios of the erosion velocity and soil permeability of the dilated zone k (i.e.

Extending XBeach to account for spatial heterogeneity of soil permeability
Alongshore (spatial) variability of dune response to surges and erosion is often controlled by three main groups of alongshore variabilities (Elsayed and Oumeraci 2016): (i) Hydrodynamic (i.e. alongshore variability of the hydraulic load in terms of water levels and wave characteristics), (ii) Topo-bathymetric (i.e. spatial variability of subaerial dune topography and fronting subaqueous bathymetry), and (iii) Geotechnical (i.e. spatial variability of soil resistance/tendency to erosion). These variability forms act individually or collectively in order to determine response of coastal barriers to storm surges and probable location of possible breaches. The mutual interaction between the hydraulic loading (L) and barrier topography (T) is satisfactorily addressed and is often described by the simulation framework of Mickey et al (2018) in Figure 3. This framework provides estimates of the mutual interaction between L and T with the XBeach model and allows examination of alongshore as well as cross-shore variabilities of storm impacts. However, such simulation framework so far omitted a specific focus on spatial variability of the soil resistance R because of incapability to quantify its order of magnitude and the incapability of XBeach to consider spatially varying soil properties within its coding framework.
Beaches, dunes, and barrier islands exhibit high levels of spatial diversity in composition and spatial heterogeneity of soil properties, including grain's D 50 , porosity, and permeability (Brain 2016;Gore et al. 2019;She et al. 2008). Accounting for such spatial variability requires mapping of each soil property individually. However, for regional studies of coastal erosion, D 50 and porosity can be assumed spatially homogeneous because the effects of their variability are significantly less important as is the effect of varying permeability (Elsayed et al. 2021;Foortse et al. 2019;Lemmens et al. 2016;Mohr et al. 2018Mohr et al. , 2021. It is hypothesized that omission of longshore permeability distribution may be one reason why most recent studies on the response of coastal barriers to storm surges (e.g. Harley et al., 2016;Muller et al., 2016;Passeri et al., 2018;van der Lugt et al., 2019) have had difficulties when comparing modeling outcomes with observations, and in determining the exact erosion volumes as well as locations, actual sizes of overwash openings and breaches.
In this work, the XBeach model is extended to account for spatially varying permeability for the first time. Varying permeability can now be supplied to the XBeach solver and considered in the analysis of coastal erosion based on Equation 5. To that end, the scalar value of the permeability in Equation 5, which was previously assumed uniform alongshore, is transformed in XBeach to a matrix of the same dimensions of the model mesh so that each cell in the mesh can be assigned a different permeability value using a new defined input file named svprfacfile.txt. This has been implemented in the XBeach code in a way similar to the spatially varying bed friction (Harter and Figlus 2017;van der Lugt et al. 2019;Passeri et al. 2018;Schambach et al. 2018). The XBeach code extensions to account for spatially varying soil permeability as well as their coding are detailed in Elsayed et al (2020). In this way, the modeling domain can be divided into subdomains and identify different permeability values for each subdomain. The above-outlined extension leads to accounting for the soil resistance to erosion R in terms of a spatially varying permeability as shown in the analysis framework in Figure 3. The interaction among the spatial variabilities of (i) hydraulic loading (L), (ii) topography (T), and (iii) soil resistance to erosion (R) can be examined more accurately. Moreover, the role of interaction among these variabilities in determining probable breaching locations can be studied. It is expected that portions of coastal dunes and barrier islands that feature low soil resistance to erosion (i.e. higher permeability) may undergo more erosion and hence breach earlier, meaning that locations of high permeability will have more potential to erosion and therefore they are more probable as breaching locations. This is further investigated by applying the extended XBeach model to the following case study.

Case study
Due to the difficulties related to the in-situ measuring of soil properties, in particular permeability, it is quite difficult to find measurements of soil properties alongshore a coastal barrier. To overcome this problem, a synthetic case study of an artificial barrier island is designed here (see Figure 4), which takes benefit of the experimental measurements of Foortse (2016). In fact, Foortse's measurements of multiple bentonite fractions and their corresponding permeability values provide the opportunity to investigate alongshore variability of resistance to erosion while knowing detailed information about the soil constituting the artificial barrier under consideration. The advantage of these measurements is that they showed only variations of the permeability with the bentonite fractions, while other soil properties, e.g. porosity, angle of internal friction and cohesion, remained almost constant (see Foortse et al. 2019). This advantage enables simulating the sand-bentonite mixtures as pure sand but with various permeability values based on the bentonite fraction. Thus, the permeability (as a representative of the bentonite fraction) remains the single parameter that influencing the soil resistance to erosion in this synthetic case study.

Description of the case study
The case study in Figure 4 represents an alongshore uniform barrier island (in terms of topography and fronting bathymetry). The alongshore uniformity of topography was intentional in order to neutralize the effect of longshore variability in topography (T) from the analysis so that the interaction between the applied hydraulic load (L) and the soil resistance to erosion (R) can only be examined. The considered domain is 400 m in the cross-shore direction and 2000 m in the longshore direction. This domain is divided from offshore to onshore into three main sections: (i) the sea section of a 300 m cross-shore width, (ii) the barrier island section of cross-shore width of 65 m that consists of a dune of a 45 m fixed base width with 20 m wide berm and (iii) the bay section of 35 m width. In the sea section, the bed level varies from (-30.00 m) at the offshore boundary to (0.00 m) at the shoreline according to the equation shown in Figure 4b.
It is assumed that this coastal zone consists of the same sand (Φ 256 µm) that Foortse (2016) used in his experiments, which has a median diameter 50 of 256 µm and 15 of 192 µm. This assumption is made to build on his detailed information about the geotechnical properties and strength characteristics of the sand, including permeability values. Five cross-shore cross-sections (Sec X1-X5) as shown in Figure 5a are set as comparison cross-sections. Besides, one longshore profile (Sec Y1) is set to pass through the dune crest of the dune. In order to study the effect of alongshore variability of permeability, it was assumed that the dune permeability is varying alongshore as shown in Figure 5a. The most left and most right 250 m (i.e. zone C (ZC) and zone C* (ZC*)) have the same value of bed permeability, i.e. 6.10E-07 m/s. The second most left and most right 500 m (i.e. ZB and ZB*) have the same value of permeability, i.e. 4.70E-06 m/s. The 500 m around the central axis of the dune (i.e. ZA) has a permeability of 4.80E-04 m/s. This means that the highest value of permeability is defined at the 500 m in the middle of the dune and the value of permeability is decreasing toward the lateral boundaries so that the variability of permeability is considered symmetric around Sec. X3. This symmetry was intentional in order to investigate the variability around the central axis (Sec X3) and to enforce a potential breach around it. Despite the variability of permeability alongshore, other soil parameters are considered spatially uniform so that the permeability remains the principal analysis parameter here.
A combination of low magnitude/high frequency and episodic hydraulic loads were applied to the barrier island synthetic model. The wind-induced high-frequency waves were defined in the model by a JONSWAP spectrum of 1.75 m significant wave height, a peak period of 6.6 s, a peak enhancement factor of 3.3, and a wave direction perpendicular to the coastline. This JONSWAP spectrum was applied along the whole simulation time (4.5 hours). The episodic load represents the storm-surge induced sealevel rise, which was considered as shown in Figure 5b. The seawater level rose from zero level (0.00 m) under the effect of both the surge and the astronomical tide and caused a surge wave. The long wave is assumed to last four hours with a peak of 6.6 m after 2 hours and varies linearly each hour.

Testing program and considered boundaries
As the considered case study have alongshore uniform topo-bathymetric properties, the interaction between the hydraulic load and the soil resistance to erosion can be considered. Alongshore uniform (long-crested) and non-uniform (short-crested) wave loads can be applied to the barrier island. Longcrested waves can be generated in XBeach by assigning a very high directional spreading parameter (s >1000). Such long-crested waves are supposed to produce alongshore uniform erosion as long as the bathymetry and geotechnical parameters are uniform alongshore. On the other side, generating directionally spreading waves produces non-uniform wave loads alongshore because of introducing wave groupings in the longshore direction. As a result, sediment exchange and dune erosion rates can vary alongshore producing different development of the foreshore in time.
Three numerical experiments (runs) were carried out, as shown in Table 1, to study the interaction between the hydraulic load and the soil resistance to erosion. To that end, the model domain was discretized every 2 m in X-and Y-directions, except at the dune location the spatial step in cross-shore direction was considered 1.0 m to enhance model performance during dune avalanching.
In Run 1, an alongshore uniform hydraulic load is applied to the barrier island after assigning a spatially uniform value of the permeability equal to 6.10E-07 m/s, which is the same value considered at ZC and ZC* in Figure 5a. The uniform permeability is assigned to XBeach as a matrix that has the same dimensions of the model mesh using the new extension of the XBeach model. Run 1 aimed at examining the interaction between alongshore uniform loading and soil properties. In Run 2, the same alongshore uniform hydraulic load in Run 1 is applied to the barrier island after assigning the spatially varying permeability values in Figure 5a. This run aimed at understanding the effect of spatially varying permeability on determining possible breaching locations and dimensions.

Soil resistance (R) represented by permeability Notes
Run 1 Uniform alongshore uniform alongshore uniform alongshore Soil permeability is considered uniform alongshore of a value equal to 6.10E-07 m/s, which is the same values at zones ZC and ZC* in Figure 5.a Run 2 Uniform alongshore uniform alongshore varying alongshore Soil permeability is considered varying alongshore the dune slices as shown in Figure 5.a Run 3 Non-uniform alongshore * The hydraulic load represents a combination of the storm surge and the JONSWAP spectrum. In Run 1 and 2, the directional spreading parameter s is assigned equal to 10000 while it is considered equal to 10 in Run 3.
In Run 3, simulation of a directionally spreading wave field of 10 degrees and spatially varying permeability were run in order to understand the interaction between spatial variability of the applied load and the spatially varying soil properties. In these three runs, the main parameters shown in Table 2 were defined to set up the model properties and its boundary conditions. For instance, the bed friction was considered uniform by setting Manning coefficient n = 0.03 [m -1/3 .s]. The uniform and spatially varying permeability are assigned to the model using the newly added keyword svprfacfile, which calls the permeability file of uniform (in case of Run1) and spatially varying values (in case of Run 2 and 3). The other XBeach parameters are kept by the defaults of XBeach and the results of the analysis are discussed below.

Results and discussion
Using the extended XBeach model, the role of uniform and spatially varying permeability was investigated under the effect of alongshore uniform and non-uniform loading by comparing the outcomes of Run 1, Run 2, and Run 3. Figure 6 shows the outcomes of these three Runs at two different times (i.e. at t = 2 and 4 h). When an alongshore uniform load is applied to alongshore uniform dune of alongshore uniform soil properties, the induced erosion was also uniform alongshore as shown in Figure 6a and b. Meanwhile, applying the same load to the same dune but with varying the permeability alongshore (as shown in Figure 5a) resulted in alongshore varying erosion (see Figure 6c and d), where locations of higher permeability values were more prone to erosion. Because of the symmetry of the permeability values in Run 2 around Sec X3 (see Fig 5a) and uniformity of the applied load, the developed erosion was also symmetric around Sec X3. One breach as shown in Figure 6d developed at the higher permeability zone (ZA), which has a permeability of 4.80E-04 m/s. The breach width equaled the alongshore width of the higher permeability zone (ZA). Less erosion developed at the zones ZB and ZB* which have less permeability (4.70E-06 m/s). Because the most left and most right 250 m (ZC and ZC*) have the same value of bed permeability of 6.10E-07 m/s in Runs 1 and 2, the predicted erosions along these zones were equal in both Runs. The relation between the permeability value and the subsequent erosion in Run 2 is shown in Figure 7a, which presents the evolution at Sec Y1 that passes alongshore through the dune crest. The dune crest was symmetrically lowering with time marching, based on the permeability value. A breach of 500 m width was developed along ZA (i.e. from y = 750 to 1250 m). The zones ZB and ZB* (i.e. from y = 250 to 750 m and from y = 1250 to 1750 m) were less eroded because they have less permeability value. The zones ZC and ZC* were less lowered because of the very low permeability. Figure 7b shows the initial and final profile (at t = 4 h) at the five cross-shore cross-sections (Sec X1-X5) and proves the symmetry of the erosion around Sec X3 when an alongshore uniform load is applied. The bed evolution at Sec X1 is the same as at Sec X5 and the same for Sec X2 and Sec X4. At Sec X3, the dune is fully overwashed, due to the higher permeability value. This means that the locations of the higher permeability are the more likely breaching locations. Figure 6e and f show the results of Run 3 when short-crested waves are interacting with the dune of spatially varying soil resistance to erosion. The alongshore varying hydraulic load applied in Run 3 induced indeed changes to the erosion pattern as compared to Run 2, in which long-crested waves of the same energy and riding on the same surge level are applied. Therefore, the erosion volumes fluctuated in the longshore direction as compared to the case of uniform loading in Run 2. Moreover, the shortcrested waves resulted in two breaches of different sizes as shown in Figure 6f. These breaches differ in size when compared to the case of uniform loading in Run 2 (i.e. Figure 6d). However, the locations of higher soil permeability (i.e. the zone ZA) are still the more prone positions to more erosion and breaching. This means that accounting for soil properties reduces the uncertainties of predicting probable locations of breaches.  Figure 8 shows a comparison between the erosion and deposition volumes induced by the uniform loading in Run 1 and 2 with the same volumes induced by the non-uniform load in Run 3. The comparison is based on dividing the study area in the cross-shore direction into 5 longshore sections: (i) from the offshore boundary to the shoreline (i.e. from 0 to 300 m); (ii) the berm section (i.e. from 300 to 320); (iii) the seaward dune front from the dune's seaward toe to the crest (i.e. from 320 to 340); (v) the bayward dune front from the dune's crest to the bay-ward toe (i.e. from 340 to 365 m); and (iv) the bay section (i.e. from 365 to 400 m).

Role of spatially varying permeability
The uniformity of the loading and soil properties in Run 1 led to uniform deposition (+ ΔV) or erosion (-ΔV) as shown by the blue lines in Figure 8. The non-uniformity of the soil properties in Run 2 led to alongshore varying erosion and deposition volumes as shown by the black line in Figure 8. The spatial variability of the hydraulic loading in Run 3 led to more fluctuations of the erosion and deposition volumes as shown by the red line in Figure 8. However, the fluctuations developed in Run 3 were still around the values achieved by the uniform loading in Run 2. This means that the variability of the soil properties had a major role in determining breaching locations, while the variability of the load is secondary. In general, offshore the shoreline represents a deposition zone due to return flow (see Figure 8a and f), where the active zone starts approximately at x = 270 m. The berm section is also a deposition zone because avalanching sand from the dune front deposit on the berm (see Figure 8b and g). The seaward dune front is an erosion zone (see Figure 8c and h). The same applies to the bay-ward dune front (see Figure 8d and i). The bay section is always a sediment deposition zone.

Soil porosity effect on alongshore variability
Dunes and other natural coastal barriers grow slowly during calm conditions as a result of aeolian and marine processes (Cohn et al. 2018). This means that such barriers are supposed to be formed of loose soils. However, in many instances, coastal barriers are so intensively human-affected and therefore either consolidated naturally during several nearshore human activities or artificially compacted to increase their strength and resiliency against coastal erosion (Brain 2016;Harris et al. 2020;She et al. 2008). Soil compaction decreases soil volume and thus increases the bulk unit weight of soil (Hanson et al. 2007;Tabrizi et al. 2017). As a result, void ratio and porosity also decrease, leading to an increase of the threshold shear stress (Mitchener and Torfs 1996;Mohr et al. 2018). As a result, permeability reduces with the increased bulk unit weight (Ranaivomanana et al. 2017). The reduction of both the porosity and the permeability hinders coastal erosion by increasing the soil resistance R (Equation 5).
When Run 2 is re-considered with decreased porosity value (i.e. 0 = 0.3 instead of 0.4), the erosion is decreased as shown in Figure 9 because of increasing the resistance and hence the modified critical Shields number according to the dilatancy theory (Equation 5). The increase of resistance arises from the fact that more densely packed sand requires more shear stress to dilate; i.e. the difference between 1 and 0 in Equation 5 increases.
As shown in Figure 9c, the effect of decreased porosity was significant at the location of higher permeability (i.e. at ZA) and vanishes at locations of very low permeability (i.e. at ZC und ZC*). Therefore, the depth of the breach, which took place at the higher permeability zone ZA, is substantially decreased with decreasing the porosity. On the sides of the breach at ZB and ZB*, the effect of decreased porosity was minor because of the lower permeability. Figure 9 depicts also that decreasing the permeability value is more efficient in mitigating coastal erosion than decreasing the porosity. In fact, porosity becomes unimportant at locations of very low permeability e.g. ZC and ZC*. Anyway, using both measures together may be more effective than using one of them in strengthening natural or manmade coastal barriers.

Summary and concluding remarks
This study addressed the alongshore variability of coastal erosion, particularly of natural coastal barriers under storm surges. As a novel contribution of this study, it introduced an extension for the XBeach model in order to account for spatial variability of permeability as a measure of spatially varying soil resistance to erosion. This extension enhanced the modeling framework of Mickey et al (2018), which is commonly used in modeling coastal morphology and mutual interactions between hydrodynamics and subsequent morphodynamics. The new extension integrates the spatial variability of soil properties into the analysis framework of coastal erosion. It, therefore, enables studying interactions among spatial variability of the hydraulic loading and induced morphodynamic while accounting for spatial variability of soil resistance to erosion. Because of the lack of field soil properties data in order to validate the model extension, this study applied the extended XBeach model to an artificial case study that took benefit from the laboratory measurements of Foortse et al (2019) for soil properties and formed an artificial barrier island of spatially varying permeability.
The analysis of the case study has revealed that soil properties of barriers, in particular permeability, control the alongshore variability of erosion to a large extent and point to possible breaching locations. The role of soil porosity was secondary, but it becomes very important at alongshore portions of the barrier that have high permeability. Moreover, low permeability portions alongshore the barrier showed very high resistance to erosion while high permeability portions resisted less and therefore breaches developed at the high permeability zones.
For purpose of applications based on this study, mapping of alongshore variability of permeability and porosity is quite important for better prediction of coastal erosion and probable breaching locations and sizes. Moreover, decreasing the values of permeability (by e.g. bentonite additives) and/or decreasing the porosity (by e.g. compaction) alongshore natural or man-made coastal defenses (e.g. levees or dikes) would function as an efficient structural mitigation measure against increased coastal erosion, breaching and subsequent significant inland discharges. Dunes, barrier islands, and sand cores of typical German dikes may be more resilient during extreme events by reducing their permeability or by increased consolidation.
The conclusions drawn in this study are based on outcomes of a synthetic case study and hence indepth experimental work that investigate the alongshore role of soil properties in determining the coastal response is a recommended topic for further research. Moreover, this study urges coastal morphologists to collect field data for soil properties to be used in future studies investigating alongshore variability of coastal erosion and examining probable coastal barriers breaching locations under extreme storm surges.

ACKNOWLEDGMENTS
The financial support of the German Research Foundation (Deutsche Forschungsgemeinschaft; DFG) for the first author (fund no. EL 1017-1/1) within the research project ReFresh (Response of Coastal Barriers and Freshwater Aquifers to Extreme Storm Surges and Flooding) is gratefully acknowledged.