Detection of Annual Coastal Boulder Mobility From 2017–2021 Remote Sensing Imagery (Gulf of Taranto, Mediterranean Sea) and Relation With Marine Weather Conditions

1. Introduction

In several coastal areas, different-sized boulders are transported and emplaced above the sea level by waves that surge over the coastline like land-crossing bores. Massive displacements of large boulders are considered morphodynamic signatures of coastal inundation events due to storms and tsunamis and can provide significant data for coastal evolution and hazard assessment [1,2,3]. From about a decade, specific bands available with commercial Earth observation satellites such as WorldView-2 allowed to define some geomorphological features of coastal deposits in detail [4,5]. The open-access to high-resolution satellite imagery via web-interface provided by Google Earth has greatly increased the possibilities of understanding the boulder formation and mobility [6,7]. Since the first availability of Landsat-8 images, online measuring tools were used to measure the horizontal dimensions of the boulders as well as their transport distance [8,9]. With the increase of the images resolution, these measurements have become increasingly reliable and are used in hydrodynamic modeling [10,11].

The boulder mobility accounts for the number of boulders displaced over a stretch of coast by water flows during a given time interval. It has attracted the attention of several scholars in recent years and has been used for various purposes among which: to test hydrodynamics equations that govern the boulder displacement [12]; to explore the ability of contemporary storms to cause geomorphological changes on rocky coasts [13]; to investigate the relationship among transport distance and storm frequency and intensity [14]. Moreover, the measure of the mobility of coastal boulders over pluriannual periods has been addressed in several studies to explore sedimentological features (cf. e.g. [15,16]).

In a previous work, by bi-temporal image analysis and geomorphological surveys, a displacement of tens of boulders occurred between July 2018 and June 2020 along the 100 km eastern coast of the Gulf of Taranto (southern Apulia region, Italy, central Mediterranean) was detected [11]. With reference to the same study area and using multi-temporal object-based analysis of 2017 to 2021 remote sensing images, in this note: additional wave-displaced boulders spread over 4 new sites are recognized; and, the annual mobility for medium to very coarse boulders was established. As a further achievement, the marine weather conditions responsible for this massive displacement phenomenon are described and discussed.

2. Study case

The Gulf of Taranto faces the north Ionian Sea for more than 330 km (Figure 1) and is characterized by a microtidal regime with an astronomical tides of 0.3 m. For such a semi-enclosed basin, maximum significant waves of 6.3 m for a 50 year return period and of 8.2 m for a 100 year return period have calculated, respectively [17,18]. Along the low lying eastern side of the gulf, jointed carbonate rocks extensively outcrop. This coast is prone to boulder formation and displacement and, because of such characteristic, two selected sites (namely Sant'Isidoro and Torre Suda, see Figure 1) are periodically surveyed with direct observation since 2017 [19,20,21]. Careful geomorphological surveys are made especially after the more severe storms.

The two 6-month climatic seasons of October–March (cold–wet, CW) and April–September (warm–dry, WD) characterize the Mediterranean climate over the Gulf of Taranto [22,23]. Indeed, the southwards migration of the Atlantic storm track and the descending mid-atmospheric throughs from Northern Europe are typical of the CW season, while good weather with scarce precipitations prevails from April to September. The Apulia region is characterized by the predominance of northerly winds, although also winds from south are significant for frequency and strength [24]. The Gulf of Taranto is exposed to strong winds generated by the seasonal prevalence of the anticyclones. Typical local storms in the central Mediterranean Sea are usually driven by middle atmosphere low pressure centers (troughs in the 500 hPa GPH weather maps) generally following a NW to SW path [23]. Their income is characterized by southerly warm currents over the south Mediterranean Sea, that can cause strong surface winds according with the position of the surface pressure minima, that can often intensify when near the middle atmosphere trough, especially in warm sea surface conditions, like in fall season. Intense southerly winds can cause quite high waves above the southern Apulia coasts, because of their potentially long geometrical fetch southwards to the Africa coasts, and in accordance with their duration [20,25].

Comparing two sets of Landsat-8 imagery, Delle Rose et al. (2021) [11] recognized 81 boulder displacements occurred between July 2018 and June 2020, and spread over 7 sites of the eastern coast of the Gulf of Taranto. One additional site (Sant'Emiliano, see Figure 1) was found along the western coast of the Strait of Otranto. Both the initial and final positions were identified for 59 displaced boulders. Instead, only the initial position was identified for 5 boulders, apparently fallen into the sea because of wave impacts. For the remaining 17 boulders, the final position was recognized while the initial one was not established because they were submerged before the displacement or were not recognized in the July 2018 image. Every remotely recognized pre- and post-displacement positions, together with the lengths of the boulder axes, were verified on the field. Meteorological investigation and hydrodynamic modeling [11], support that both the storm Vaia (28-30 October 2018, see [26] for the synoptic framework) and the storm Detlef (11-13 November 2019, see [27] for the synoptic framework) can be considered as the main causes of the massive boulder displacement. However, it cannot be ruled out that some boulders may have been displaced under rough sea conditions due to favorable combinations of hydrodynamic processes (cf. e.g. [28,29]).

Aerial images taken during July 2019 at an average flight altitude of 180 m, have recently made available on the site "Territorial Planning Service" [30] (Servizio Assetto del Territorio, in Italian), which is managed by the governmental Authority of Apulia. In this paper the eastern coast of the Gulf of Taranto is visually examined by a computer workstation using four sets of Landsat 8 images (dated July 2017, July 2018, June 2020, and September 2021) and the aforementioned July 2019 aerial images to detect the boulder mobility at the annual scale. The visual interpretation was performed according to the rudiments of Geographic Object-Based Image Analysis (GEOBIA) [31,32,33], and selected images were elaborated in an open-source vector graphics editor.

3. Method

3.1. Multi-Temporal (Geographic Object-Based) Analysis

The interpretation of remote sensing imagery is a visual problem-solving activity which explores the relationship between objects present in the image to obtain the recognition of the image content. It is first based on the knowledge and experience of the analyst [34,35,36]. In geosciences, like aerial photo-interpretation, visual analysis of satellite imagery relies on several qualitative features of Earth's surface image, that are tone, texture, pattern, shape, context, and scale [37,38]. Reflectance properties of the surface materials, lighting conditions, and image processing and display, determine the tone of the objects, namely their colour and brightness. Tonal differences between the objects is usually the first observed characteristic in image interpretation. The texture, as a combination of the magnitude and frequency of tonal change, is related to the scale and resolution of the images. The spatial arrangement of tones and textures which make up the image scene results in the pattern, and is mainly related to geomorphology, topography, vegetation, and human imprint. The image characteristics change depending on the time of image acquisition, due to changes in lighting conditions, vegetation cover and soil moisture content [39].

The generation of geographic information from remote sensing imagery that enable users to perceive, interpret and respond to specific issues is the main purpose of GEOBIA. This analytical perspective first aims to develop theory, methods, and tools for each issue, to finally replicate and exceed experienced human interpretation of the images in semi-automated or even automated manner [31]. A geographic object is defined as a "bounded geographic region that can be identified for a period of time as the referent of a geographic term" and, as such, can be represented in the map [40]. The temporal bound well accounts for mobile objects such as displaced boulders during the image acquisition time. The property which endows coastal boulder with an enduring identity, and which differ it from the surroundings is the detachment from the parent rock. However, the demarcation of objects like the boulders in remote sensing images can be rather problematic because coastal areas are structured into continuously varying patterns which sometimes hamper the boundary drawing at a given scale (cf. e.g. [41]). Again, the accuracy of the boulder demarcation depends on the degree of illumination and the shadow effects within the images, which are results of the sun-object-sensor geometry.

In the visual interpretation of satellite imagery aimed to the detection of coastal boulders, special attention is given to the masking effects of vegetation cover and topographic shadow. They can abruptly reduce the chances of recognizing objects in images, especially the shadow of a cliff cast on a coastal plain [7]. On the other hand, shadows have long since played an important role in remote sensing, the effects of shadowing having been used to determine ground features and object height in applied sciences and archaeology [42,43]. Nevertheless, brightness and textural characteristics of satellite images are exploited in geosciences mainly to detect tectonic lineament, soil feature, and landslide contour [44,45,46]. Finally, it must be noted that multi-temporal satellite images analysis by Google Earth Engine has been recently tested in improving landslide visibility [47].

Boulder usually is a complex solid [48,49,50]. However, in the first stage of the geomorphological investigation, it is considered an idealized-shaped rectangular cuboid, and the size described by the dimension of its axes, a (major), b (middle), and c (minor). The axes dimensions can be used to establish the Flatness Index (FI) = (a+b)/2c and the Shape (Sh) of the boulders [51,52]. According to the Udden-Wentworth grain-size scale, boulders range from the class of fine boulders, with a middle axis >0.25 m, to the class of medium blocks, with a middle axis of up to 16.4 m [52]. The middle axis ranges from 1 to 2 m for coarse boulders, and from 2 to 4.1 m for very coarse boulders. In remote sensing analysis, the two axes of the boulder up-face presented in the images can be measured. They are designated as the apparent major (a_r) and medium axes (b_r), with the minor axis being usually perpendicular to the image [8,9,10]. Actually, since the boulders may be inclined rather than horizontal, a_r and b_r are equal or smaller than the major and middle axis, respectively. Dealing with boulder displacement, other geometrical features that can be determined by the computer workstation are the initial (x_i) and final distances (x_f) of the boulder from the coastline, and the related Transport Distance (TD). In this study, TD was determined by estimating the change in horizontal position of the center point of the boulder in the images. With reference to this latter, the geographical coordinates of each boulder position were carried out.

To verify the reliability of boulder detection from remote sensing imagery, geological field surveys were made. The on-site investigations aimed also to examine the features of the ground covered by boulders before the displacement. A boulder displaced from a given initial position leaves a temporary trace on the ground surface, more or less similar to its planar shape. Where the boulder is detached from the parent rock, the trace is named socket and consists in a fresh and not weathered rock surfaces that is initially not covered by lichen. The socket may still be detectable several years after the displacement [53]. Instead, the trace left by a boulder that has already been uprooted from its geological substrate and has been lying on the ground for some time, usually has a shorter duration and is detectable in satellite images under favorable sunlight conditions [11].

Following the approach of previous studies and using Google Earth Engine [7,10,54], to detect the displacement of boulders with apparent a-axis of at least 1 m, visual comparisons were made among four sets of Landsat-8 images at an eye elevation of 50 to 200 m. The sets of the used images are dated: 9 July 2017, 20 July 2018, 28 June 2020, and September 2021 (the day of month is not reported by the resource), respectively. The object (boulder) demarcation was manually performed taking care to consider the shadowing conditions. As a matter of fact, the shape of the boulders can apparently change according to the shadows cast at different times of the day [55]. Because of resolution limits of the aerial images analysed by the regional viewer service [30], the geometrical characteristics of the boulders were usually less clear than in the Landsat-8 imagery. However, the quality of the aerial images was sufficient to compare shape and position of the displaced boulders with those shown in the satellite images. Actually, similar techniques have already been used successfully in previous studies (cf. e.g. [55,56]). Given the above, in this study the 2017-2021 annual mobility of medium, coarse, and very coarse boulders spread over the study coast was detected.

3.2. Evaluation of the Marine Weather Conditions

To describe the marine weather conditions in the eastern Gulf of Taranto during the 2017-2021, data extracted by publicly available data sets have been used. Considering the purpose of the study, special attention was given to the storminess conditions. A sea storm is ''a succession of sea states in which the significant wave height exceeds a fixed threshold for a duration of at least equal to 12 h'' [57]. Moreover, to be relevant by the geomorphological point of view, a storm must have ''the potential to significantly alter the underlying morphology and expose the backshore to waves, currents and/or inundation'' [58].

A screening of the Central Mediterranean 2017-2021 storms has been made using the web archives of the GLOBO-BOLAM-MOLOCH model cascade [59]. To identify the middle-atmosphere pressure conditions that allow to check the presence of stormy conditions, the 500 hPa geopotential height maps of the GLOBO model archive were careful screened. Then, in order to obtain information about the offshore wind field (fetch length and direction, average wind speed at 10 m height, storm duration) the archives of the mesoscale models BOLAM (10 km horizontal resolution) and MOLOCH (3 km horizontal resolution) were used. The storms were selected considering a minimum duration of 12 hours together with a minimum fetch usually greater than 500 km, besides a minimum wind speed of 12 m/s. These are typical conditions achieved by the storms affecting the Gulf of Taranto, and correspond to significant wave heights generally above 3 m (sea state ≥ 5 of the Douglas Sea Scale).

An evaluation of the characteristic wave height H₀ and the wave period T in the sea storms has been made starting from the wind characteristics (duration R, length of the fetch over the sea surface F, fetch direction D, and wind speed U) as established by the forecast map archive of the BOLAM-MOLOCH model. The wind characteristics were used to calculate the spectral peak wave height from the similarity relation equations [60,61]. They relate T, H₀ and R to the wind speed at 10 m height offshore U and to the fetch over the sea F:

$$gT/U=0.286{(gF/U^2)}^{1/3},$$

(1)

$$g{H}_0/U^2=0.0016{\left(gF/U^2\right)}^{1/2},$$

(2)

$$gR/U=68.8{\left(gF/U^2\right)}^{2/3},$$

(3)

where g is the gravity acceleration [61] (p. 117).

The duration R of stability of wind conditions, estimated from the wind maps, has been used as a limiting factor to calculate an effective fetch F for the equations 1 and 2, whenever the actual geometric fetch appears to be longer. It must be noted that H₀ is nearly equivalent to the significant wave height [62,63].

Data of wind gauge stations placed at Santa Maria di Leuca, Gallipoli, and Porto Cesareo (Figure 1 for location) are below used to infer the nearshore wind conditions during storms presumed to be responsible for changes in boulder position (Section 5). Due to the incompleteness and inhomogeneity of the data sets, these conditions can be only partially defined. The Santa Maria di Leuca and Gallipoli stations are managed by the Apulia Civil Protection that publishes comprehensive annual reports (named Annali Idrologici, Hydrological Annals) [64,65,66,67,68]. These report the daily average and maximum together with four 30 min averages at four selected hours of the day (00, 06, 12, 18). However, the Gallipoli station has been in operation only since November 2020. The Porto Cesareo station belongs the Apulia Meteomarine Network (Sistema Informativo Meteo Oceanografico delle coste Pugliesi, in Italian); its data are online publicly available and consist of 15 min measures of wind speed and direction [69]. The Porto Cesareo station has been in operation since February 2019.

An estimation of the nearshore wave energy related to the boulder displacements is here attempted. The theoretical framework is briefly described in what follows. The displacement of boulders due to wave impact on the coast starts when the hydrodynamics forces exceed the resistance forces of net friction [70,71]. Nandasena et al. (2011,2013) [1,72] and Nandasena (2020) [73] elaborated a set of equations to calculate the minimum flow velocity (V) required to set in motion and onshore transport the boulders, depending on Pre Transport Setting (PTS) and Movement Type (MT). The PTS refers to some physical conditions (scenario) of the boulder before the displacement. The boulder can be isolated (sub-aerial, SA, or submerged, SB), laterally limited by joints (joint-bounded, JB), located at the cliff-edge (CE) [1,72]. These scenarios determine the possible three MTs (saltation/lifting, sliding, and overturning/rolling) during the boulder displacement. SA and SB boulders can be lifted, slid, or rolled by the water flow; due to the movement constraints, JB boulders can only be lifted, while CE boulders can be lifted or rolled [70,73].

The magnitude of V can give an estimation of the energy of nearshore waves of the storm responsible for the change in position of the boulders [56,74,75]. The incipient motion formulas are given next [1,72,73]; for SA and SB boulders moved by saltation/lifting, sliding, and rolling, the minimum flow velocity is respectively:

$$V\ge {\left[\frac{2({\mathrm{p}}_s/{\mathrm{p}}_w-1)gc\cos \mathrm{\theta }}{{\mathrm{C}}_L}\right]}^{1/2},$$

(4)

$$V\ge {\left[\frac{2({\mathrm{p}}_s/{\mathrm{p}}_w-1\left)gc\right(\mathrm{\mu }\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{\theta }+\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{\theta }}{{\mathrm{C}}_D\left(c/b\right)+\mathrm{\mu }{\mathrm{C}}_L}\right]}^{1/2},$$

(5)

$$V\ge {\left[\frac{2({\mathrm{p}}_s/{\mathrm{p}}_w-1\left)gc\right(\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{\theta }+\left(c/b\right)\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{\theta })}{{\mathrm{C}}_D\left(c^2/b^2\right)+{\mathrm{C}}_L}\right]}^{1/2},$$

(6)

For JB boulders moved by saltation/lifting, the minimum flow velocity is:

$$V\ge {\left[\frac{2({\mathrm{p}}_s/{\mathrm{p}}_w-1\left)gc\right(\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{\theta }+\mathrm{\mu }\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{\theta })}{{\mathrm{C}}_L}\right]}^{1/2},$$

(7)

For CE boulders moved by saltation/lifting and rolling, the minimum flow velocity is respectively:

$$V\ge {\left[\frac{2({\mathrm{p}}_s/{\mathrm{p}}_w-1)gc}{{\mathrm{C}}_L-\mu {\mathrm{C}}_D\left(c/b\right)}\right]}^{1/2},$$

(8)

$$V\ge {\left[\frac{2({\mathrm{p}}_s/{\mathrm{p}}_w-1)gc}{{\mathrm{C}}_L-{\mathrm{C}}_D\left(c^2/b^2\right)}\right]}^{1/2},$$

(9)

In the Eqs. 4, 5, 6, 7, 8, and 9: b = b-axis length (m); c = c-axis length (m); ρ_s = density of boulder (kg m^-3); ρ_w = density of seawater (kg m^-3); g = acceleration due to gravity (9.81 ms^-2); θ = bed slope angle (degrees); μ = coefficient of static friction along ground surface; C_L = lift coefficient; C_D = drag coefficient.

One can observe that c is present in all the equations, thus, where hydrodynamic information are to be obtained from displaced boulders, the c-axis length must be measured on-site. On the other hand, using the geographical coordinates obtained from Google Earth, the studied objects can be easily found in the field with a global positioning system (GPS). Alternatively, using specific bands available with commercial Earth observation satellites, c-axis length can be successfully retrieved from remote data [5]. In the course of the field survey, the inclination to the horizontal of the boulders I and the dimensions of their axes (a_f, major, b_f, middle, c_f, minor) have to be measured. Moreover, the Lithology (Li) of the boulders (which, in turn, allows to determine ρ_s) and clues on PTS and MT can be obtained by field observations [20,53]. The data collected by satellite imagery and geological field survey for the present study are shown in Section 4 and Appendix A.

For the Eqs. 4, 5, 6, 7, 8, and 9, the values of the coefficients C_D and C_L are crucial in determining the reliability of the results [76]. However, for a given site, it is difficult to establish the most appropriate value of C_D among those experimentally calculated [70,77], especially because the influence of environmental factors like the substrate roughness and flow turbulence. Nevertheless, there are no specific investigations performed to calculate C_L for the boulders [78], while hydrodynamic site-specific measurements should even be necessary to correctly reach the target [72]. Despite such issues, the values of 1.95 for C_D, and 0.178 for C_L are commonly used in literature [1,11,56,75,79]. About the μ coefficient, the value of 0.7 is the most commonly used. In this paper, in order to make comparisons with previous studies possible, the latter values are used (see Section 4).

4. Results

The 2017 to 2021 annual boulder mobility for the eastern coast of the Gulf of Taranto and the related marine weather conditions are reported below. In particular, the new recognized displacements are described in Section 4.1}, while the multi-temporal check performed on the displacements recognized by [11] is shown in Section 4.2. The main features of the more severe storms and the estimation of the nearshore wave energy are given in Section 4.3.

4.1. New Detection of Boulder Displacement

Multi-temporal analysis allowed the recognition of 38 boulder displacements occurred during the investigated time period and spread over 4 previously undetected sites (in red in Figure 1). In Table 1, these sites are reported from north to south. For 25 out of 38 boulders, both the pre- and post-transport positions were recognized, thus the TD was established.

The results obtained from images analysis are shown below, while numerical data and supplementary field measures and observations are in Appendix A. The initial and final coordinates of the displaced boulders are in Tables A1,A3,A5,A7. The lengths of the boulder axes measured on satellite imagery and on-site, respectively, together with other features (initial and final distances from the coastline; transport distance; lithology; shape; flatness index; pre transport setting) of the displaced boulders are in Tables A2,A4,A6,A8.

The northernmost site, Santa Caterina, is quite far from the other three sites (about 25 kilometers), while these latter are close together (a few kilometers from each other, Figure 1). Moreover, the stretch of coast of Santa Caterina is characterized by a nearly vertical cliff about 5 m high unlike the others which do not present abrupt slope breaks.

The 2 boulders displaced at Santa Caterina (Tables 1,A1,A2) experienced different movements (Figure 2). SCa was positioned at the cliff-edge at least up to July 2019. In its place, since June 2020, a socket is present (Figure A1 in Appendix A) while the detached boulder is not detectable in the images. Likely, it fell back onto the sea bottom after the impact of a high-energy wave. It must be noted that the socket is hardly visible in September 2021 image, probably because unfavorable sunlight condition. The boulder SCb is clearly distinguishable respect the surrounding in July 2017 image, differently from the subsequent image. SCb was still in the same position in July 2019, while it has a new position in both June 2020 and September 2021 images. Currently, it lies on a shrub vegetation area (Figure A2 in Appendix A). The displacement of SCb is quantifiable in a TD of about 15 m (Table A1) from SW to NE (Figure 2f). Both the SCa and SCb boulders were displaced between July 2019 and June 2020 (Figure 2c,d).

Along the stretch of coast of Pilella (Figure 1 for location), 14 displaced boulders were identified by multi-temporal analysis. For 11 of these boulders both the initial and final positions were detected and, thus, the TD calculated (Table 1). Middle axis ranges from 1 to 1.5 m, except for one case (PILj) which is nearly 2 m (Table A4); TD ranges from 1.5 to 5 m (Table A3) with vectors oriented about SSW-NNE. The initial positions of the boulders apparently do not change in the July 2017, July 2018 and July 2019 images; likewise, the final positions do not change in the June 2020 and September 2021 images. However, it was not possible to determine whether PILc boulder changed position before or after July 2019. By way of example, the northernmost boulders of Pilella site are shown in Figure 3. Comparing the boulder shadow casts of July 2017 and July 2018 images, detailed planar shapes of PILl, PILm, and PILn can be inferred. Given this accurate localization, the boulders can be easily identified in the low quality July 2019 image. Two boulders of this site are inclined of 55–60° with respect to the ground (Table A4). Such a geomorphological setting can make axis measurement from remote sensing images rather unreliable. As an example, for the PILn boulder (see Figure A3 for its arrangement on the coast), b_r resulted a third smaller than b_f (Table A4).

In June 2020 and September 2021 images, PILl, PILm, and PILn are overturned in comparison with the positions of the previous images. Their exposed surfaces, not covered by lichen, are whitish and easily distinguishable from the surrounding gray texture. Consequently, they had been overturned by water flow, roughly oriented SSW-NNE, between July 2019 and June 2020.

The site of Posto Rosso, where 16 displaced boulders were identified (Table 1) is located about 1.5 km southeast of Pilella (Table A5). For 6 of these boulders the initial position was not established, while the final position was established for all the 16 boulders. Their middle axis ranges from 1 to 2 m (Table A6), while the calculated transport distance from 2 to 17 m (Table A5); vectors are oriented about SSW-NNE. In Figure 4, the analysed time series remote sensing images of the central sector of the stretch the Posto Rosso coast are reported. As in the previous site, due to the different sun-boulder-sensor geometries of July 2017 and July 2018 images, the planar shape of the boulders can be defined in detail and, consequently, recognized also in the low quality July 2019 image. The pre-displacement positions of ROe, ROg, ROh, and ROl boulders have not been identified. Where they were in SA PTS, the characteristics of the analyzed images did not allow their detection. Alternatively, these boulders were in SB PTS before to be displaced. Except for ROe boulder, the changes in position occurred between July 2019 and June 2020 (Figure 4c,d), while no changes apparently occurred between July 2017 and July 2018 (Figure 4a,b), and July 2018 and July 2019 (Figure 4b,c); ROl changed position between June 2020 and September 2021 (Figure 4d,e; see also Figure A4b in Appendix A for the geomorphological setting). The change in position of ROo and ROp boulders occurred between July 2018 and June 2020 as detected by the satellite images. Unfortunately, these boulders are hardly recognizable in the 2019 aerial image, thus the time interval during which they moved cannot be further detailed.

Among the 4 new sites identified, Posto Rosso is the one with the greatest number of displaced boulders, which are large in size. These calcarenite boulders were detached from a surface that corresponds to the stratigraphic contact with the limestone (see Figure A4a, A4c,d in Appendix A). Apparently, the stratigraphy is here crucial for the boulder formation.

The southernmost site identified by multi-temporal analysis is Torre San Giovanni (Figure 1 for location). For 3 out the 6 boulders displaced, the pre-displacement position and TD have been determined (Table A7). Their middle axis ranges from about 1 to more than 1.5 m (Table A8). The transport vectors are up to almost 10 m in length and approximately S-N oriented, as shown in Figure 5f. All the boulders of this site experienced position changes between July 2019 and June 2020 (Figure 5a,b,c,d,e). The SGc and SGd boulders have apparently finished the movement because stopped by other boulders lying on the coastal platform (Figure A5 in Appendix A).

As a whole, visual analysis of the 2017-2021 time series remote sensing images allowed to confidently establish that at least 33 out 38 boulders of the new detected sites have been displaced between July 2019 and June 2020. For 25 of them the TD was determined (Table A9). One boulder (SGd) changed position between July 2018 and July 2019; another boulder (ROl) between June 2020 and September 2021. The remaining three (PILc, ROo, and ROp) are hardly identifiable in the 2019 aerial image; anyway, they experienced the change in position between July 2018 and June 2020 as established by the comparison of the satellite images.

4.2. Check on Previously Detected Displacement

Given the above finding, the multi-temporal investigation was performed on the 81 displaced boulders previously detected by bi-temporal images analysis [11], with the aim to detect their annual mobility for the 4 years considered herein (see Sections 1,2). In what follows, some features of the boulder displacement are given with reference to the sites of Punta Prosciutto, Torre Suda, Mancaversa, and Punta Pizzo (see Figure 1).

Eighteen boulders that have changed position between July 2018 and June 2020 were previously recognized along Punta Prosciutto headland [11]. Actually, by multi-temporal analysis have been found that all the 18 boulders have changed position between July 2019 and June 2020. For 13 boulders both the initial and final positions were established, thus the TD was determined (Table A10). Several changes in positions are shown in Figure 6. The case of the cluster composed by the PRi, PRj, PRk, and PRl boulders is particularly enlightening. Before the displacement, they were some meters to the southwest and differently arranged in comparison with the post-displacement arrangement (Figure 6f). One SSW-NNE transport vector of about 5 m is drown in the figure; actually, each of the boulders has been transported from about 4 to 6 m. The boulder PRk is now vertically positioned, thus it shows middle and minor axes in remote sensing images. The other 3 boulders of the cluster lie inclined on the coastal plateau.

The boulders PRq and PRr appear the result of the breaking of a singular boulder that likely occurred during the transport (Figure 6a,b,c). Again, this geomorphological process took place between July 2019 and July 2020 (Figure 6c,d). Moreover, PRr was overturned between July 2020 and September 2021, thus changing again position (Figure 6d,e).

The largest displaced boulder detected by Delle Rose et al. [11] was found at Torre Suda (Figure 1). This boulder (named SUa) sizes 5.4 x 4.6 x 1.9 m and has been displaced by about 9 m from SSE to NNW (Figure 7). It must be reported that, due to its middle axis dimension, SUa belonging to the class of the fine blocks according to the modified Udden-Wentworth classification (cf. [52]). It apparently slid and rotated about 180° over a submarine terrace plain just below the sea level, constrained in moving from the low cliff from which it detached (Figure 7f). The constraint explains the anomalous trajectory. SUa is clearly identifiable in the aerial image taken in July 2019 (Figure 7c) where it is placed in the same position of the July 2017 and July 2018 satellite images (Figure 7a,b), while its socket is not distinguishable both in June 2020 and September 2021 images because it is close to the sea level.

By multi-temporal analysis has been found that 11 out the 14 displacements detected at Torre Suda [11] have evidently occurred between July 2019 and June 2020, although one of the displaced boulders (the SUi boulder, see Appendix A) had already undergone a first position change between July 2018 and July 2019. For the remaining three boulders (SUj, SUk, and SUl) was not established whether they were moved before or after July 2019 because of their small transports (less than 1 m).

Apart from SUa, other 9 boulders detected by [11] have middle axes that exceed 2 m in length. Due to such a feature, they belonging to the class of the very coarse boulders [52]. Among these, the case of the MAa boulder is instructive for the present study. This boulder was detected at the Mancaversa site (Figure 1) and sizes 3.4 x 2.1 x 0.5 m. Due to an unfavorable sun-object-sensor geometry, MAa is few distinguishable in the July 2018 image (Figure 8b). As a matter of fact, its presence can only be inferred from the narrow shadow cast which marks its NE side. Differently, in both the July 2017 and July 2019 images, the planar shape of MAa is well defined by topographic shadow. This boulder was displaced between July 2019 and June 2020 (Figure 8c,d), apparently sliding and rotating about 150° over the coastal plain (Figure 8f). Differently from the case of the SUa boulder (see above), the movement of MAa was only partially constrained by topographic features. After new field investigation (Appendix A) MAa is herein considered a SA boulder, differently from the previous evaluation (cf. [11]).

Differently by the cluster composed by the PRi, PRj, PRk, and PRl boulders (see above), other clusters were found with few changes in mutual spatial arrangement of the boulders after the displacement. The case of cluster composed by the PIa, PIb, PIc, PId, and PIe boulders are considered below (Figure 9). This cluster appears evenly moved by about 5 m from SSW to NNE (Figure 9f). Apparently, no collision between the boulders occurred.

Sixteen out the 18 displacements found at Punta Pizzo [11] have occurred between July 2019 and June 2020 as detected by multi-temporal images. For two boulders (PIm and PIo) was not established whether they were moved before or after July 2019.

In conclusion, by the check on the displacements detected by [11], it is found that no less than 76 out of 81 boulders have been displaced between July 2019 and June 2020. Moreover, one of these 76 boulders had already undergone a first position change between July 2018 and July 2019, while another has still changed position between June 2020 and September 2021. The remaining 5 boulders have not been identified in 2019 image due to their small size and/or transport distance (less than 1.5 m), both inconsistent with the low resolution of the aerial image. However, as for 3 boulders reported in Section 4.1 they experienced the change in position between July 2018 and June 2020.

4.3. Storms and Boulder Incipient Motion

The main features of the identified storms affecting the Gulf of Taranto during the considered 4 years (see Section 3.2) are reported in Table 2. All the selected storms occurred during the CW seasons. Marine weather conditions with sea state 5 happened in 4 days of the 2017-2018 CW season. The 2018-2019 CW season has been characterized by only two severe storms (including Vaia) both with sea state 6. Four storms with sea state from 5 to 7 (the highest value belonging to Detlef) hit the study area during the 2019-2020 CW season. Finally, 2 storms both with sea state 6 occurred during the 2020-2021 CW season (Table 2).

As mentioned above, an estimation of the energy of the nearshore waves responsible for boulder displacement can be obtain by a careful use of the hydrodynamic equations 4, 5, 6, 7, 8, and 9. However, these equations are affected by some issues (see Section 3.2). As regards the values of the lift, drag, and static friction coefficients, the ones nearly constantly adopted in coastal boulder deposit studies (0.178, 1.95, and 0.7, respectively; cf. [76]) are here used to allow comparisons with previous studies. The values used for the density coefficients are: 1.8 for calcarenite boulders and 2.2 for limestone boulders; 1.02 for the seawater. Data on a dozen of the biggest boulders (7 coarse and 5 very coarse boulders) displaced between July 2019 and September 2020 are used for the assessment (Table 3).

5. Discussion and Conclusions

5.1. Annual Mobility of Coastal Boulders

Along the eastern coast of the Gulf of Taranto, one hundred and nineteen boulders that changed position between July 2017 and September 2021 have been detected by remote sensing image analysis (Section 4.1, 4.2). Moreover, two boulders (SUi and PRr) were apparently moved twice during this lapse of time. No boulder was displaced between July 2017 and July 2018. One boulder (SUi) was certainly displaced between July 2018 and July 2019. At least 110 displacements (including the second transport of the SUi boulder) occurred between July 2019 and June 2020 (Table 4). The PRr boulder experienced its second displacement between June 2020 and September 2021, that is the same interval of the change in position of ROl boulder. For 8 boulders was not possible to establish if they experienced the transport between July 2018 and July 2019 or between July 2019 and June 2020 due to the low resolution of the 2019 aerial image (Sections 2, 4).

Even if "coastal boulders may be mobile under a range of wave conditions" [29], the transport and emplacement of a number of coarse and very coarse boulders is the morphodynamic signature of inundation and erosion processes due to tsunami events or exceptional wind-wave conditions [6,70,73]. During July 2017- September 2021, a western Mediterranean earthquake caused very small waves in the Gulf of Taranto unable to determine boulder displacements [11,80]. Consequently, the cause of the massive boulder displacement occurred during the interval 3 (Table 4) must be found among the storms occurred during the 2019-2020 CW season (Table 2). Nevertheless, with reference to previous statement [11], storm Vaia have to be ruled out as possible causative storm.

5.2. Causal Storm Inferences

Four severe storms hit the eastern side of the Gulf of Taranto during the 2019-2020 CW season (Tables 2, 4). The 12-13 November storm attracted the interest of several weather bureaus. It was named "Detlef" by the Meteorological Institute of Berlin Free University [81], and such name was adopted in some studies [3,82,83] (in [20] storm Detlef is provisionally named "Ionics19"). However, it was also named "Trudy" and "Bernardo" by the US National Aeronautics and Space Administration, and the Spanish Meteorological Agency, respectively [84,85,86]. On 10 November, a wide trough elongated from the Northern Europe toward the Western Mediterranean creating a center of low atmosphere pressure. The cyclone transformed into a subtropical storm by 0600 UTC the 11 November [87] and then moved eastward over the Central Mediterranean. Later, it reached the North Adriatic Sea causing an exceptional sea level height in the Venice lagoon [23,88,89]. Regarding the effects of storm Detlef on the study area, it must be noted that a strong wind over the eastern side of the Gulf of Taranto started soon after midnight of 12 November from the southeast, and persisted for 24 h from almost the same direction with an offshore fetch of about 900 km and an average wind speed of 20-22 m/s. The calculated characteristic wave height is 8.7 m (Table 2), that is even greater than the maximum significant wave for a 100 year return period (Section 2).

A similar synoptic condition caused the 24 November storm, however, with a weaker depression in the middle atmosphere. Again, the strongest winds blew from the south-eastern direction over the southern Apulia coast, with an intensity of 18-20 m/s (Table 2). In this case, the fetch was of about 800 km, while the duration of these conditions of about 12 h. Both the 22 December 2019 and the 2 March 2020 storms were caused by two middle atmosphere lows coming from the north-eastern direction with respect to the Italian peninsula [11]. This difference caused a different impact over the study area, with more variability of the wind speed and direction. In the case of December 2019, the prevailing direction was from the western direction over the southern Apulia coast, with a shorter fetch because of the geographic shape of the Gulf of Taranto (Figure 1). The BOLAM model shows a wind of about 14-16 m/s at 10 m over the sea, associated to a duration of 24 h and a fetch of about 400 km. In the case of March 2020, the strongest wind blew from the South, with a lower intensity of about 12-14 m/s, a duration of 12 h and a fetch of about 600 km. The calculated characteristic wave height for the storms occurred between the 24 November 2019 and the 2 March 2020 ranges from 3.4 to 4.9 m (Table 2), values lower than that of the maximum significant wave for a 50 year return period (Section 2).

Searching for the main cause of the massive displacement phenomenon occurred during the third interval time of Table 4, complementary items came from the stretches of coast that are surveyed with on-site geological observations and measures (Section 1). By surveys made after storm Detlef, 6 and 11 displaced boulders were found at Torre Suda (SUi, second displacement, SUj, SUk, SUl, SUm, SUn) and Sant'Isidoro (SIc, SIf, SIg, and 8 additional boulders that are not detected by remote sensing imagery analysis because their small size or short displacement), respectively [19,20,21] (Figure 1 for sites location). It must be noted that after storm Vaia, only one case of boulder displacement at Torre Suda site were detected (SUi, first displacement, see Figure A6a,b in Appendix A), while no displacements were found at Sant'Isidoro site.

The strong power of storm Detlef on the eastern coast of the Gulf of Taranto was also shown by the caused damage [20]. By using the Storm Power Index (SPI) introduced by Dolan and Davis (1992) [90], the intensity of the storms reported in Table 2 can be calculated. SPI is defined as the "storm’s duration times the square of maximum significant wave height" [90]. A number of studies have shown a strong correlation between SPI and the wave impact on the coasts (see e.g. [91,92,93]). Considering that H₀ is nearly equivalent to the significant wave height [62,63], the storm intensities can be assessed as reported in Figure 10.

Considering the obtained values of SPI, the quite different ability of the two most severe storms to determine morphodynamic signatures in the study area is surprising. It can be guessed that the interplay between wind speed, and length and duration in time of the fetch involved, can determine a sort of threshold-like behavior between storms that do not cause massive changes (including Vaia) and the most intense one (Detlef). Moreover, by the timely distribution of the boulder displacements (Table 4), it is inferable that the storm of 22 December 2019 may have caused the change of position of some boulders. On the other hand, the 3 December 2020 storm can be considered as the most probable cause for the two displacements occurred during the interval 4. Anyway, advanced hydrodynamic modeling based on non-linear equations to also consider the TD (Tables A2, A4, A6, A8, A9, A10), could better shed light on this physical process [94,95,96].

5.3. Nearshore Wind Conditions and Wave Energy

The nearshore wind conditions during the storms of the interval 3 (Table 4), can be inferred from 2 out 3 of the wind gauge stations placed along the eastern coast of the Gulf of Taranto (Figure 1). Unfortunately, Gallipoli station was not operational during the 2019-2020 CW season (Section 3.2). The data carried out from the datasets of Santa Maria di Leuca and Porto Cesario stations are reported in Appendix B. They confirm the stronger wind conditions occurred during storm Detlef in comparison with the ones occurred during the other storms of the 2019-2020 CW season. The 12-13 November 2019, the wind speed 30 min averages measured at Santa Maria di Leuca station ranged between 15 and 18 m/s for about 24 hours, with a peak value of 24 m/s. Wind speed and duration during other storms were lower (Table B1). The measurements of the Porto Cesareo station confirm the above. The wind speed was around 20 m/s for almost the first half of November 13, while only during December 22 did it exceed 15 m/s for about 4 hours (Figure B1, B2, B3, and B4).

Usually, wind speed is significantly higher over the sea than land while wind direction is affected by more or less apparent changes in crossing the coastline. Land station measurements of wind speed can be halved in comparison with offshore values [97,98]. However, it must be noted that the Santa Maria di Leuca station has an elevation of 26 m a.s.l. and is placed about 100 m inland, while the Porto Cesareo station has an elevation of 12 m a.s.l. and is placed few tens of m inland. Since the two stations are close to the coastline, the anemometric values reported above can be considered as indicative for the nearshore conditions. The consistency with the wind speed values extracted from the GLOBO-BOLAM-MOLOCH model cascade (Section 4.3, Table 2) confirms the above.

The minimum flow velocities required to set in motion the boulders can be carefully used to assess the nearshore wave energy (Section 3.2). Regarding the study case, for the biggest boulders displaced during the 2019-2020 CW season, V ranges from 2.7 to 8.1 m/s (Table 3). Exceptional coastal storms should produce flow velocity at shore below 5 m/s [74,99]. Apparently, the displacement of some boulders (e.g., ROn, SGe) required an unusual amount of energy for flows produced by wind waves. Values of V up to 14 m/s were calculated to explain boulder displacements due to inundation events driven by typhoons in the Asia-Pacific region [74,79,100,101,102]. Instead, for the Mediterranean region, V values up to 12 m/s are reported in literature to displace SA or SB boulders [1,56]. The results of this study therefore appear to be placed in a coherent framework. However, for the use of the data in coastal hazard assessment, an improvement of the method is suggested below.

Equations (4), (5), (6), (7), (8), and (9), could provide some unreliable minimum flow velocities to set in motion the boulders due to different causes [50,76]. The first concerns the values to be attributed to the coefficients C_D and C_L. The drag coefficient can be affected by environmental factors and actually ranges from 1.2 to 5.0 [70,76]. On the other hand, as established by calculation and numerical simulations [78,103], the lift coefficient can be about an order of magnitude greater than the commonly used value (that is 0.178, see Section 3.2). Such an issue needs to be adequately faced as the research continues. Another set of problems concerns the real shape of the boulders and the lengths of the axes. The boulders are usually complex solids, while using the measured lengths in the equations (see Section 3 for the method; Tables A1, A3, A5, and A7 for the values), they are implicitly considered as rectangular cuboid. Consequently, the flow velocity could be overestimated because the difference between theoretical and actual size. To overcome such an issue a new approach that consider the real dimensional features of the displaced boulders, obtained by modern geophysical survey techniques, has been recently proposed [50].