Uncertainties in shoreline position analysis: the role of run-up and tide in a gentle slope beach

Abstract. In recent decades in the Mediterranean Sea, high anthropic pressure from increasing economic and touristic development has affected several coastal areas. Today the erosion phenomena threaten human activities and existing structures, and interdisciplinary studies are needed to better understand actual coastal dynamics. Beach evolution analysis can be conducted using GIS methodologies, such as the well-known Digital Shoreline Analysis System (DSAS), in which error assessment based on shoreline positioning plays a significant role. In this study, a new approach is proposed to estimate the positioning errors due to tide and wave run-up influence. To improve the assessment of the wave run-up uncertainty, a spectral numerical model was used to propagate waves from deep to intermediate water and a Boussinesq-type model for intermediate water up to the swash zone. Tide effects on the uncertainty of shoreline position were evaluated using data collected by a nearby tide gauge. The proposed methodology was applied to an unprotected, dissipative Sicilian beach far from harbors and subjected to intense human activities over the last 20 years. The results show wave run-up and tide errors ranging from 0.12 to 4.5 m and from 1.20 to 1.39 m, respectively.


Introduction
Mediterranean beaches are well known for their high environmental, economic and sociocultural value. In the last few decades, most of these beaches have been subjected to demographic growth from increasing tourism and commercial activities (Cooper et al., 2009). To support these activities, new defence structures have been built along some beaches, and although these structures have reduced local erosive effects, they have also increased erosion on neighboring coasts (e.g., Griggs, 2005;Stancheva et al., 2011;Manno et al., 2016). Coastal erosion is a relevant problem that involves both socioeconomic resources and private properties, and its assessment has long been an issue of international interest involving political decision-makers and researchers (Douglas and Crowell, 2000;Phillips and Jones, 2006;Anfuso et al., 2011;Rangel Buitrago and Anfuso, 2015). Historical beach evolution, erosion, and the retreat/accretion of shorelines have been analyzed using aerial and satellite images (e.g., Thieler et al., 2009;Fletcher et al., 2003;Genz et al., 2007;Anfuso et al., 2011;Dolan et al., 1980Dolan et al., , 1991. Each remote image is often used to represent a year, and therefore the "shoreline" position identified and digitalized from each image becomes representative of all shoreline positions in that specific year. The Coastal Engineering Manual (U.S. Army, 2008) defines "shoreline" as the intersection between land and water body, but to choose a suitable proxy that accounts for the spatial and time variability (Bush et al., 1999), this boundary must be localized. Among different shoreline proxies (Boak and Turner, 2005), the wet-dry boundary is clearly identified in aerial images by the different colors of sand during the drying process. Because it is more sensitive to run-up fluctuations than astronomical tide variations (Dolan et al., 1980), the wet-dry boundary is a stable shoreline proxy that has been applied by several authors for various applications regarding localization and analysis of shorelines (e.g., Pajak and Leatherman, 2002;Moore, 2000;Moore et al., 2006;Stockdon et al., 2002;Robertson et al., 2004). Thieler et al. (2009) developed a method to assess the beach evolution trend by means of aerial imageries, implemented in a soft-ware extension to ESRI ArcGIS© v.9+, the Digital Shoreline Analysis System (DSAS), that can calculate the shoreline rate-of-change statistics starting from multiple historical shoreline positions. This method has the advantage of considering uncertainties due to positioning and measurement errors . The positioning errors are strictly connected to physical phenomena and can affect the analysis precision because the erroneous position of a shoreline is assumed to be "actual" for the considered year. Uncertainties from tides and wave storms (seasonal variability) are linked to the "exact position" of the shoreline during the aerial shooting , whereas measurement uncertainties are linked to errors of image processing and digitizing conducted by technicians who identify and map the shoreline position for several observation years . Several authors (Genz et al., 2007;Rooney et al., 2003;Romine and Fletcher, 2012) used DSAS to evaluate both positioning errors and measurement errors, neglecting the error due to wave run-up and astronomic tide fluctuations. By contrast, other authors (e.g., Virdis et al., 2012;Manca et al., 2013) added to the positioning uncertainty the effects due to wave motion, calculating the run-up by means of the empirical formula of Hunt (1959). In this paper, an interdisciplinary method that more accurately assesses shoreline positioning error caused by wave run-up and tidal fluctuations in DSAS analysis is presented. Wave run-up was calculated using a numerical model cascade, which includes a wave spectral model and a shallow water propagation model. Tide effects were evaluated using the daily variation of astronomic and meteorological tide. With this method, a dissipative sandy beach of the western coast of Sicily (Italy) was analyzed, an interesting case study because, in the last decades, it has been heavily impacted by human activities. This beach represents a practical case in which accurate identification of the shoreline position with extreme fluctuations is fundamental to forecasting inundation areas or planning effective beach management practices.

Methodology
The methodological goal was to better evaluate positioning errors caused by wave run-up and tide for DSAS applications, an ArcGIS extension used to compute the shoreline rate of change (Thieler et al., 2009). The latter was evaluated by five different methods to compare the related results. The first method considered the end point rate (EPR), calculated by dividing the shoreline shift by the time elapsed between the oldest and most recent shoreline position. The second method used the linear regression rate (LRR) of change based on the determination of least-squares regression lines of all the shoreline points of each transect. The third method used a weighted linear regression (WLR), in which the weight w is a function of the variance of the mea-surement uncertainty (Thieler et al., 2009): where e is the shoreline uncertainty value. The fourth and fifth methods are based on the analysis of distances rather than rates. The fourth method considers the "net shoreline movement" (NSM), the distance between the oldest and youngest shoreline positions for each transect, and the fifth considers the "shoreline change envelope" (SCE), the distance between the farthest and closest shorelines to the baseline at each transect. To assess the total uncertainty (σ T ) affecting each shoreline position, the following relationship was assumed (Virdis et al., 2012): where the uncertainty σ i is the standard deviation of the itype error; σ d is the digitizing error determined by digitizing several times the same feature on the image; σ r is the orthorectification error, considered as the root mean square error (RMSE) for photogrammetric blocks; σ co is the image coregistration error arising from the RMSE of misalignment between single pixels from the set of images obtained by the rectification; σ p is the pixel error assumed equal to the pixel size; and σ wr and σ td are, respectively, the wave run-up and the tide errors estimated in this study (discussed later). Note that the first four errors are related to intrinsic characteristics of the used images, how they were taken and how they were processed, whereas the last two are related to specific geomorphologic, mareographic, and wave characteristics of the beach examined. Variables σ td and σ wr represent position errors that may result in noticeably higher values than the others; therefore, special care is required during their evaluation, which is the focus of this study.
To improve the evaluation of the wave run-up and tide uncertainty (σ wr , σ td ) with respect to the use of empirical formulas found in the technical literature (e.g., Virdis et al., 2012;Manca et al., 2013), various mathematical models were applied. A hydraulic study, conducted on the basis of a geomorphologic study, determined the effects of wave motion and tide fluctuation on the shoreline position. To this aim, offshore wave parameters were used to simulate wave propagation from deep water to run-up on the beach, whereas a tidegauge dataset was used for analysis of tide fluctuation. The whole mathematical process for the run-up calculation can be summarized by the following steps: (a) select offshore buoy dataset collection; (b) propagate waves from deep to intermediate water by means of a wave spectral model; (c) generate random waves from a JONSWAP spectrum; (d) propagate waves from intermediate water up to the swash zone with a Boussinesq-type model; and (e) conduct run-up analysis. This mathematical process has been validated using in-field measurements as described in Sect. 4. The Boussinesq-type model considered in the present paper is able to propagate the waves from a relatively shallow water depth (kh = 0.7 where k is the wave number and h is the local water depth) up to the shoreline. Assuming a JONSWAP spectrum, the significant wave height and wave period were converted into the time series of an energetically equivalent irregular wave train, which was then propagated using the shoreline Lagrangian numerical model of Lo Re et al. (2012b). In the shoreline model, a Boussinesq-type model for breaking waves with the governing equations solved in the ζ − u form was implemented, where ζ is the free surface elevation and u is the depthaveraged horizontal velocity. The values of the variables ζ and u were calculated inside the wet domain, whereas the shoreline position (defined by means of its horizontal coordinate ξ(t) perpendicular to the coast) and its velocity u s were calculated by means of the Lagrangian shoreline equations. In the case of an orthogonal wave attack such as the one considered here, the variable ξ is only a function of time, i.e., ξ = ξ(t), and the kinematic condition at the shoreline is the following: Such a relation states that the fluid particles at the shoreline remain along the shoreline. Moreover, the momentum equation at the shoreline that must be also be considered in order to close the problem in dimensional form reads where ∂ζ /∂x| s is the derivative of the surface elevation evaluated at the shoreline and F fric is the bottom friction force, evaluated as follows: in which h is the local depth and f is the bottom friction coefficient. When the value of F fric becomes too large, due to the small value of the total water depth, a threshold is used. In such a case, the dependency on the water depth is eliminated and the bottom friction is assumed to be only a quadratic function of the depth-averaged velocity: where C f is a coefficient that was assumed equal to 5.0 m −1 . The propagation of the offshore wave characteristics to shallow water was carried out by the well-known SWAN spectral propagation model Holthuijsen et al., 1993;Ris et al., 1999). The SWAN results obtained for the 5 m bathymetric line were then used as input for the Boussinesq-type model by Lo Re et al. (2012b) which, coupled with a specific Lagrangian model for shoreline movement, allowed simulation of wave swash and run-up. The wave run-up error σ wr was finally estimated by analyzing the resulting shoreline movement over time. Note that the offshore wave parameters were the only source in the SWAN propagation model. Moreover, because the SWAN simulation covers a small region in intermediate water, the wind input, the wave drop due to white-capping, and the wave drop due to bottom friction were not considered (Rusu, 2009). Indeed, in large deep water regions with very shallow water, the wind and the bottom friction could play a significant role and cannot be neglected. For wave propagation by SWAN, a 2-D unstructured grid following the Delaunay rule was implemented, constructed in accordance with Monteforte et al. (2015) using a density function in which the triangle sizes depended on local water depth and wavelength. The node elevation was calculated by a linear interpolation of bathymetric data from nautical charts. The Lagrangian model used for shoreline movement discriminates between wet and dry regions to simulate run-up and run-down along surveyed sections. For the ith section the standard deviation of the horizontal shoreline movement over time was calculated by where S wr, i is the standard deviation of the vertical shoreline fluctuation computed by the model, and tan α i is the section slope. Finally, the wave run-up error for the whole beach σ wr was estimated by σ wr = σ 2 wr, 1 + σ 2 wr, 2 + . . . + σ 2 wr, n (n − 1) .
The tide uncertainty σ td was assessed by processing the tide measurements recorded by a mareographic station. For each year with measures, the standard deviation of the tide measurements, S td , was first computed, and the standard deviations of the horizontal tidal fluctuations of the same sections used for run-up assessment were then evaluated using an equation formally identical to Eq. (7). The tide uncertainty for the whole beach, σ td , was finally assessed with the same equation used for the run-up error (Eq. 8).

The case study: Marsala beach
The case study of the dissipative beach ( Fig. 1), Lido Signorino, extends in a north-south direction for about 3.5 km between Cape Torre Tunna and Cape Torre Sibilliana. Its slope ranges between 1.5 and 10.8 • and the direction of beach exposure (Fig. 1) is about 140 • , between northwest and south-southwest. The Egadian Islands, in particular Favignana, shield the beach in the 320 • N direction. The geographical fetch is limited from the west by the Spanish coast, from the south by the African coast, and from the northwest by the Sardinian coast.
The buoy belongs to the Italian Wave Buoy Network (RON) and the rose, obtained by processing available data  The beach suffers from intense anthropogenic use, especially houses placed too close to the shoreline (Fig. 2), which has caused progressive destruction of dunes and their associated natural supply. In the first 50 years of the 20th century, the dunes were uniform from north to south and about 5 m high, whereas today they are discontinuous, about 2.5 m high and mainly located in the less developed southern area.

Validation of run-up assessment by means of field measurements: a Marsala beach
In order to validate the whole mathematical process used for run-up assessment, field measurements were performed (Lo Re et al., 2012a). The wave run-up on sandy beaches can be measured in several ways depending on the general aim and on the amount of detail required.
Records of the shoreline positions can be obtained by resistance run-up gauges or by video cameras. The technique adopted in the present paper is based on a high-frequency video monitoring system (Holman and Sallenger, 1985). Such a kind of technique allows the acquisition of several images by means of a digital video camera. The choice of the position of the camera was a fundamental task because the camera has to shoot the whole studied area but at a short distance, in order to obtain the maximum level of detail from the recorded images. In particular, positions of the swash were measured on a transect across the beach, normal to the shore (Fig. 3). For this transect a line was built using a rod at 0.5 m intervals. The first stake was a piezometer and it was next to the beach berm. The second stake of the line was placed at a distance of 5 m from the piezometer. The line stakes on the beach profile were georeferenced using control points from a previous topographic survey. The video camera was placed at a distance of 10 m from the line of stakes (orthogonally), and it was used to record 240 min continuously. The shot videos were digitized in order to extract the wave run-up of each wave. When a wave reached a stake, the data were recorded. The horizontal run-up distance was calculated starting from SWL obtained from the water level inside the piezometer. Finally the corresponding run-up value was estimated by considering the beach profile. Each run-up measurement (R) was recorded in time windows of 30 min (eight windows in total) according to Nielsen and Hanslow (1991). For all recorded data, the Rayleigh distribution was fitted by using the leastsquares method. The application of the Rayleigh distribution to our data allowed estimation of the 2 % run-up (R 2 % ).  The expression of the Rayleigh cumulative distribution function is reported as follows: in which R 100 is the value transgressed by 100 % of the waves, i.e., the lower limit of the distribution, and L zwm is the vertical scale of the distribution, i.e., the shape parameter. Moreover, to perform such a validation, wave parameters from the buoy of Mazara del Vallo were used (Wave Buoy Network, managed by the Italian Institute for Environmental Protection and Research, ISPRA). In particular: (1) significant wave heights, H s (m); (2) peak period T p (s); and (3) mean wave direction D m ( • N). The extraction time period goes from 11:30 to 15:30 of 29 March 2011. The waves shown in Table 1 correspond to the sea states recorded by the buoy half-hourly.
The obtained wave run-ups are reported in Table 1. Such a table also shows, for each time window, the results obtained with the empirical formula by Nielsen and Hanslow (1991).
The R 2 % run-up determined by means of the Rayleigh distribution of field measurements is also shown.
The analysis of the 2 % run-up (R 2 % ) highlights that both methods give acceptable results. In particular, the numerical model has an average percentage error of 2.81 % and the empirical formula gives an average percentage error of 18.11 %. The numerical Boussinesq model gives overall results closer to the field measurements and, for this reason, it was chosen for simulations of wave run-up in this study.

Method application and results
Five orthorectified aerial images were used to assess time variations of the shoreline position during the 1994-2007 time span (Table 2). Each image was georeferenced (WGS84 -UTM 33N) by 6-10 evenly spaced ground points. For each observation year, these images were used to form a photomosaic covering the whole coast studied. The shoreline relating to a photo-mosaic was traced and digitized manually using the wet-dry proxy, as suggested by Virdis et al. (2012) for dissipative beaches of the Mediterranean Sea. In order to reconstruct wave conditions at the day the aerial photos were taken (Table 4), data recorded from the Mazara del Vallo buoy were analyzed (Fig. 4). In particular, 3-hourly wave parameters were used for a total of 40 sea states for 5 days (each day has its specific beach profile, wave and tide). Every single sea state (see Table 4) was then propagated throughout the numerical domain by SWAN in stationary mode (Fig. 4). At the 5 m bathymetric line, the wave spectrum output of SWAN was used to generate a wave time series (Table 4), which in turn was used as input to the Boussinesq model to assess the wave run-up in 26 sections. Table 3 shows the average slope for each day analyzed. Therefore, 1040 (26 sections × 40 sea states × 5 days) near-shore simulations were conducted for each offshore wave.
The considered sections (Fig. 4) are distinguished from the transects (discussed later) by an S preceding the number. In addition to the 26 sections, 68 transects orthogonal to the present shoreline (Fig. 4) were generated (with the DSAS application) at about 50 m from one another to better analyze the shoreline changes and the related erosion/accretion rates between the transects themselves.
The Boussinesq wave propagation on 26 sections produced a number of run-up values that were then processed, obtaining the 2 % wave run-up (R 2 % ) and the run-up standard deviation σ wr, i (i = 1, 2, . . ., 26) relating to each survey day (Fig. 5). Comparison between the run-up values and the standard deviation indicates, for each year, consistency be-tween the R 2 % and σ wr, i trends (Fig. 5). The wave run-up errors σ wr (Eq. 8) were summarized in Table 6, together with the other uncertainties.
For each day with measurements, the standard deviation of the tide measurements (S td ) was first computed and then the standard deviations of the horizontal tidal fluctuations of the 26 sections were evaluated by Eq. (8) ( Table 5). As is well known, tide fluctuation measurements include the meteorological effects.
Based on the uncertainties (Table 6), the five shoreline rate-of-change indexes mentioned earlier were evaluated using DSAS for each of the 68 transects. Indexes WLR, EPR and LRR (Fig. 6) and indexes NSM and SCE (Fig. 7) were plotted, with positive index values indicating shoreline accretion and negative values indicating recession.

Comparison of shoreline rates of change
The trend differences in the five indexes showed that variations among WLR, EPR, and LRR (Fig. 6a) are consistent, exhibiting generally similar trends on both accretion and recession. Several accretion zones are clearly distinguishable in transects 12-13, 33-40, 47, and 65-67. The accretion rate is about 0.6 m year −1 for transect 47, 0.51 m year −1 for transects 12-13, and 1.0 m year −1 in transects 33-40, 47, and 65-67. By contrast, a recession is evident for transects 56-64, with a rate about −2.5 m year −1 , as well as for transects 48-50, at −1.6 m year −1 . For transects 20-22 and 40-43, lower recession rates of −1.18 and −1.28 m year −1 , respectively, are observed. In contrast to the other indexes, NSM and SCE often present opposing trends (Fig. 6b), in which a relative maximum of 1 may correspond to a relative minimum of the other. This pattern depends on the different definitions of the two indexes, but disparate conclusions can be drawn if one or the other criterion is adopted. Indeed, SCE represents the total change in shoreline movement for all available shoreline positions and is not related to their dates. Figure 6b shows that the NSM trend on the whole is consistent with the other indexes' trend ( Fig. 6a), except for transects 36-37 and 44-45. The lowest shoreline index values (stable areas) were generally localized in transects from 3 to 9 (8 local maximum), from 22 to 34 (24 local maximum) and from 50 to 55.     In particular, in transects 12-13 a local deposition cusp was detected. The studied beach is not characterized by a rhythmic morphology, but massive presence of beach cusps could produce erroneous results as mentioned by Anfuso et al. (2016).
Transects 51 and 55 proved to be affected by even lower shoreline movements, despite considerable anthropogenic disturbances of the dunes, this for all the indexes shown in Fig. 6a, b. Unlike the SCE index, the NSM trend (Fig. 6b) showed a reliable average equilibrium between beach accretion and recession, although in transects 56-64 a noticeable general recession occurs, in accordance with EPR, WLR, and LRR results (Fig. 6a).
Finally, the NSM indexes relating to each period between two consecutive surveys were compared to highlight the shoreline evolution at each transect along the whole period 1994-2007 (Fig. 7), showing that accretion periods alternated with recession periods at each transect. In a few transects (e.g., transects 11 and 19) NSM was effectively unchanged, whereas in many others it changed noticeably from one period to another. Moreover, during 2006-2007, accretion prevailed over recession along the entire beach, whereas recession was prevalent during 2000-2006, and during 1994-1999 and 1999-2000, accretion and recession basically balanced each other. The different behaviors (higher or lower accretion or recession) of different beach stretches observed in a given period followed the specific beach conformation of the stretches as well as the presence of Posidonia oceanica leaf deposits. Note that all the indexes considered detected higher shoreline changes in transects 58-68, where vast deposits of Posidonia oceanica leaves were present.

Conclusions
To analyze beaches by aerial images utilizing DSAS or similar applications, technicians usually neglect positioning uncertainties or assess them by means of empirical formulas without the use of hydraulic models. In this study, a new approach was adopted that assesses positioning uncertainties by analyzing wave motion and tide effects.
The hydraulic models used in our methodological approach consist of a nearshore model (SWAN) and a more accurate dispersive model (Boussinesq) that provides a more accurate description of the hydrodynamics in the nearshore area, a fundamental step to estimate the oscillation of the shoreline. The proposed method can be used in all types of coastal areas (steep beaches, gentle slope beaches, etc.) and can reproduce the hydrodynamics of a large area, not just the hydrodynamics of one point. Moreover, the models can reproduce almost all wave propagation effects (diffraction, refraction, reflection, shoaling, breaking, etc.). The diachronic analysis of Lido Signorino shows a low shoreline variability, except for the southernmost coastal stretch that has a high variability due to anthropic and natural causes.
The methodology adopted here provides high accuracy in wave run-up calculation, resulting in a more accurate σ wr error estimation as highlighted from comparison of the present model with in situ run-up measurements (see Sect. 4).
Using an application like DSAS that neglects or underestimates the σ wr error may prevent determining whether a beach is in erosion or accretion, and a retreat rate close to the total uncertainty would not be acceptable. Furthermore, method accuracy is valuable in beach monitoring and management, especially when more sustainable methods are needed to sustain coastal resources. An integrated management of the coasts must be interdisciplinary and consider the dynamic process of the beaches, mainly when the beaches are largely urbanized and anthropized.