# Assessing and modeling the dynamics and persistence of mussels in rocky-shore microhabitats

Romina Vanessa Barbosa, Cédric Bacher, Fred Jean, Marion Jaud, Jérôme Ammann, Yoann Thomas

► To cite this version:

Romina Vanessa Barbosa, Cédric Bacher, Fred Jean, Marion Jaud, Jérôme Ammann, et al. Assessing and modeling the dynamics and persistence of mussels in rocky-shore microhabitats. Marine Environmental Research, 2025, 210, pp.107321. 10.1016/j.marenvres.2025.107321. hal-05168067

HAL Id: hal-05168067

https://hal.univ-brest.fr/hal-05168067v1

Submitted on 18 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Assessing and modeling the dynamics and persistence of mussels in rocky-shore microhabitats

Romina Vanessa Barbosa<sup>a,b</sup> <sup> DID </sup>, Cédric Bacher<sup>c,*</sup> <sup> D </sup>, Fred Jean<sup>a</sup> <sup> D </sup>, Marion Jaud<sup>d</sup> <sup> D </sup>, Jérôme Ammann<sup>d</sup> <sup> D </sup>, Yoann Thomas<sup>a</sup> <sup> D </sup>

## ARTICLE INFO

### Keywords:

Species distribution

Intertidal height

Disturbance

Resilience

Remote sensing

Habitat mapping

Photogrammetry

Markov chain model

## ABSTRACT

In intertidal rocky shores, sessile organisms are subject to local variations due to a complex interplay of physical and biological drivers. Here, we examined the patterns and the dynamics of the distribution of a key benthic ecosystem-engineer species, the mussel *Mytilus galloprovincialis*. We characterized areas of loss, gain, and persistence of mussel space occupancy using drone images with a resolution of 1 cm across two rocky shores spanning 3121 m² and 3499 m². Then, we assessed the effect of microhabitat slope, microhabitat orientation, and intertidal height on the characteristics of mussel distribution for the two rocky shores. Notably, substantial losses ranged from 56 % to 80 % depending on the rocky shore. The hierarchical impact of habitat topographic features on mussel occupancy diminishes in the order of rocky shore, intertidal height, and microhabitat features. To unravel the dynamics further, we calculated persistence, resilience, and the time to recover from disturbance using a Markov Chain model. Persistence time decreased at the extremes of the rocky shore —both its top and bottom edges. Contrary to expectations, mussels at these extremities exhibit higher resilience, preventing their disappearance despite significant losses. Our results provide insights into the structural dynamics and emphasize the need to consider long-term responses to environmental changes. This study underscores the potential of integrating models with drone technology for monitoring intertidal populations, offering a powerful tool to comprehend and anticipate the enduring impacts of environmental shifts on spatial distribution.

## 1. Introduction

Understanding and predicting changes or the persistence of species distribution are major goals in population ecology related to important issues (Silliman et al., 2011; Wootton and Forester, 2013). In marine ecosystems, the persistence and changes in the distribution of species are generally addressed at high spatial scales where major drivers, such as those related to climate change, modify species occurrence and geographical boundaries (Poloczanska et al., 2013; Thomas et al., 2020; Curd et al., 2023). However, the interaction processes between organisms and their environment take place first at the individual level (i.e., ecophysiological processes) and thus on a small scale. This likely has repercussions on species distribution dynamics that may contradict large-scale trends due to local conditions that do not necessarily follow large-scale trends (Helmuth, 2002).

In intertidal rocky shores, sessile organisms are subject to local variations due to a complex interplay of physical and biological drivers, such as microclimate temperature and small-scale topographic features, and interspecific competition (Dayton, 1971; Seabra et al., 2011; Meager et al., 2011; Choi et al., 2019). For instance, Meager et al. (2011) showed that topographic features and temperature influence the abundance and body-size distribution of invertebrates on a rocky shore. The relationship between environmental conditions and the spatial distribution of species seems evident, although the particular role of each factor and the driving mechanisms is complex and not well understood (e.g., Barbosa et al., 2021). In these habitats, mussels are sentinel species studied to understand the effect of environmental conditions on species' dynamics and distribution (Helmuth et al., 2011, 2016). Mussels are one of the dominant species inhabiting the mid-intertidal zone of temperate rocky shores worldwide (Paine, 1974; Paine and Levin, 1981) and are

<footnote><sup>a</sup> Univ Brest, CNRS, IRD, Ifremer, LEMAR, IUEM, Plouzané, F-29280, France</footnote>

<footnote><sup>b</sup> University of Victoria, BC, Canada</footnote>

<footnote><sup>c</sup> Ifremer, DYNECO, Plouzané, F-29280, France</footnote>

<footnote><sup>d</sup> Univ Brest, CNRS, Ifremer, Geo-Ocean, IUEM, Plouzané, F-29280, France</footnote>

<footnote>* Corresponding author. Ifremer, DYNECO, Plouzane, F-29280, France.</footnote>

<footnote>E-mail address: cedric.bacher@ifremer.fr (C. Bacher).</footnote>
considered paramount ecosystem engineers (Borthagaray and Carranza, 2007; Arribas et al., 2014). The dynamics and persistence of these ecosystem engineer species are of main interest because their presence increases biodiversity (Arribas et al., 2014; Scrosati et al., 2021), and, thus, changes in their distribution can produce cascading effects on community structure and productivity (Silliman et al., 2011).

Merging field observations and models allows the simplification of complex systems by integrating the multiple processes involved in the dynamics of species distribution. Such integration provides the tools to understand spatial and temporal changes, as well as to test hypotheses (Guichard et al., 2001; Wootton, 2001). The development of high-resolution spatial analysis from aerial images has allowed the mapping of intertidal species distribution with high accuracy and detail (Guichard et al., 2000; Gomes et al., 2018; Garza, 2019; Barbosa et al., 2022), and could help to understand its long-term dynamics. The majority of studies on mussel distribution dynamics focus on observations of recovery after experimental mortality, such as in-situ removal of mussels in a specific square area (e.g., Carrington et al., 2009; Guichard et al., 2003). Hunt and Scheibling (2001) studied natural mortality by disturbance events on two rocky shores, but their conclusions were based on the entire area rather than on the microhabitat scale. Studying the distribution dynamics using drone images would not only allow to map entire rocky shores but also to determine the spatial patterns among microhabitats. Combining observations and models has been successfully applied to assess the colonization process or successions between species (Wootton, 2001; Guichard et al., 2003). Markov chain approach is based on probabilities of transitions between different states, e.g., the probability of transition from an occupied to an empty state and vice versa. Among other advantages, Markov chains allow the assessment of the equilibrium state and estimate the time needed to reach this equilibrium after a disturbance, and offer the prospect of a change in scale (Hill et al., 2002). Experimentally tractable local scales can be scaled up to recreate larger-scale patterns (Wootton, 2001; Guichard et al., 2003).

In a previous study, drone images were used to map mussel distribution (Barbosa et al., 2022). Results showed that the distribution of mussels was primarily related to the intertidal height gradient and, to a lesser extent, to habitat morphology in slope and orientation (Barbosa et al., 2022, see also Gilek et al., 2001). Based on that pattern, we hypothesized that the persistence of mussel aggregations determines higher occupancy at intermediate intertidal heights. According to this hypothesis, the occupation rate (proportion of occupied area relative to the available area) would increase if mussels persist over time. In contrast, the top and bottom intertidal boundaries should present a higher probability of loss with a consequent decrease in occupancy. Similarly, we hypothesized that the topographic microhabitat conditions, slope, and orientation, are related, but to a lesser extent, to mussel dynamics.

Here, we used both mapping (drone surveys) and modeling (Markov chain modeling) to assess the effect of microhabitat features on the dynamics of the distribution of a benthic intertidal mussel, *Mytilus galloprovincialis*, and to estimate key properties of the temporal dynamics of mussel distribution. We studied two rocky shores of the same bay as a natural laboratory where the general oceanographic conditions affect microhabitats similarly within both shores but differ depending on the shore orientations (mesoscale). Specifically, we evaluated the effect of shore and microhabitat intertidal height, slope, and orientation angle on the mid-term (annual changes) dynamics of *M. galloprovincialis* distribution. We also assessed the long-term effect of intertidal height on the persistence and resilience of the population, as well as the potential role of mortality measured as loss area on the species' occupancy pattern along the intertidal height gradient. To do that, we 1) evaluated the temporal dynamics of the distribution (loss, persistence, and gain) of mussel aggregations during a period of a year from drone images, 2) related the observed dynamics with the microhabitat intertidal height, slope, and orientation angle, and 3) modeled the long-term dynamics of loss and gain areas to estimate their long-term persistence, resilience

capacity, and time to return to equilibrium after cumulative mortality through the year.

## 2. Methods

### 2.1. Study site

The study was conducted at the “Le Petit Minou’ site, located on the western coast of French Brittany, France (Fig. 1A). This site is characterized by a macrotidal regime at the entrance of the Bay of Brest. The site comprises two rocky shores on both sides of a sandy beach, referred to as the West and East shores based on their positions. The West shore is oriented to the south-southeast. In contrast, the East shore is oriented to the west (Fig. 1B). The West and East rocky shores encompass a total available intertidal area within the mussel distribution range (3.32–5.80 m) of 3121 m² and 3499 m², respectively (Barbosa et al., 2022).

The “Le Petit Minou’ site is characterized by a dominance of *Mytilus galloprovincialis* (Simon et al., 2020), although hybrid individuals are common in the Brittany region (Bierne et al., 2003). The mussel bed structure is mainly formed by monolayer aggregations on both rocky shores (Barbosa et al., 2021).

### 2.2. Drone image acquisition and processing

Drone surveys were performed on two dates, June 5, 2019, and June 23, 2020, on each of the two rocky shores. Images were taken at 30 m altitude with a DJI® Phantom 4 pro (DJI Sciences and Technologies Ltd., Shenzhen, China) with autonomous flight mode during low tide. Images were acquired with an 80 % overlap to optimize further image analysis. Several ground control points (GCPs) were included in the surveys and consisted of red discs distributed in the study area. The position of each GCP was recorded with a centimeter precision using a real-time kinematic (RTK) GPS. The environmental conditions of light and sun position were similar on both dates, which resulted in similar contrast image quality (Fig. S1). Images were processed to construct the digital surface model (DSM) and orthomosaics of both rocky shores and dates using PhotoScan v.1.4.0 software (Agisoft, 2018). The GCPs were tagged in the corresponding images manually and used during the image processing to build the orthomosaic and the SDM. We employed the Time-SIFT approach (Feurer and Vinatier, 2018), which involves aligning all images from different periods (specifically, June 2019 and June 2020) and then separating the previously aligned images for each date to create 3D dense point clouds and generate the SDM and orthomosaic. This process ensured perfect alignment between the orthomosaics of both dates, facilitating a more precise evaluation of temporal changes (Feurer and Vinatier, 2018; Feurer et al., 2024). Due to the high number of cells in the resultant layer (cell resolution 0.5 × 0.5 cm) and the difficulty of processing the total ~16M cells/rows, we resampled the layer at a cell resolution of 1 × 1 cm in QGIS v.3.14.16 (QGIS Development Team, 2020) by averaging neighbor pixels for further analyses. Subsequently, we identified changes at each cell of the orthomosaic, i.e., at a centimeter resolution.

### 2.3. Identification of mussels' distribution on two dates

The identification of mussels was carried out by calculating the Mussel Visualization Index (MVI) from the orthomosaics (Barbosa et al., 2022).

$$MVI = \left( B - \frac{(R + G)}{2} \right) 100$$

where B represents blue, R represents red and G represents green, the three bands (RGB) of each cell (pixel) in the orthomosaics. This index increases the differentiation of blue mussels from the surrounding areas
<figure><img src="image_3.png" /><figcaption>Fig. 1. A) Study site location, “Le Petit Minou”, France, and B) the West and East rocky shore areas surveyed by the drone in “Le Petit Minou”.</figcaption></figure>

by enhancing their characteristic blue color. The MVI index showed good accuracy for differentiating mussels from empty rocky areas under mostly uniform light conditions, i.e., without big shadow areas, especially when mussel cells are not mixed with algae, like in our study site (Barbosa et al., 2022). Once the MVI was calculated, we applied a date-and site-specific threshold to differentiate mussels from other species and substrates based on user-selected sampled pixels of mussels and non-mussel areas from the constructed orthomosaic (details in Barbosa et al., 2022).

Areas with the presence of algae and lichens in the bottom and upper intertidal height (outside the mussel range) as well as some shadow areas exhibited similar MVI values to mussel areas. To enhance mussel identification, we performed the following workflow: 1) mask by hill-shade distribution, 2) identify mussel intertidal height range, 3) select site's cells within mussels' MVI range, and 4) select mussels' cells by intertidal height range (after Barbosa et al., 2022). This workflow was applied for each date, with the dataset and analyses from June 2019 identical to those presented in Barbosa et al. (2022). The intertidal height range of mussel distribution in the study site was recorded as between 3.2 and 5.80 m above the Lowest Astronomical Tide (LAT) in 2019 (Barbosa et al., 2022). Cells outside this range were considered false positives (false mussel coverage cells) and discarded from both orthomosaics (2019 and 2020).

The accuracy of the resultant distribution was assessed by sampling 130 cells from each orthomosaic (total of 520 cells; 130 cells × 2 rocky shores × 2 dates). For each orthomosaic, 30 cells corresponded to mussels and 100 to non-mussel cells. Non-mussel cells were manually selected from the orthomosaic, where high resolution allowed unequivocal visual differentiation of mussels from the substrate or other species' cells. Mussels' cells were also manually selected to represent the wide spatial distribution of mussels. Mussels and non-mussels' reference cells were compared with the result of the mussel identification process in a confusion matrix. The confusion matrix presents the proportion of cells corresponding to mussels and non-mussels' cells that were correctly identified, as well as the proportion that was erroneously classified. We reported the classification accuracy of mussels and non-mussels' cells separately. Then, the overall accuracy (Congalton, 1991) and the Kappa index (Kruskal and Wallis, 1952; Cohen, 1960), were calculated to determine the overall accuracy of the supervised classification. Overall accuracy represented the percentage of cells correctly classified, whereas the Kappa index represented the concordance between the supervised classification and the real distribution of mussels. The Kappa index can take values from -1 to 1, with -1 equivalent to zero concordance, whereas 1 indicates perfect concordance (true positives and true negatives).

## 2.4. Characterization of the dynamics of mussel's occupied area

The observed temporal dynamics of mussel distribution were derived by comparing both orthomosaics, from June 2019 and June 2020, to create a layer with three categories: 'persistence', gain' and 'loss'. Areas initially covered in 2019 that continued to be occupied in 2020 were classified as 'persistence areas'. Areas initially covered in 2019 that showed an absence of mussels in 2020 were classified as 'loss areas', and areas with an initial absence of mussels that were occupied in 2020 were classified as 'gain areas'. To characterize the distribution dynamics, the average and standard deviation (SD) of net area of persistence, gain, and loss on each rocky shore was measured, and represented as proportions relative to the entire available area (3.2–5.8 m above the LAT) on the shore. Additionally, the observed probability of loss, gain, and persistence was calculated along the intertidal height range on both rocky shores. Loss probability was defined as the ratio of occupied cells in 2019 that became empty in 2020; Persistence probability as the ratio of occupied cells in 2019 that remained occupied in 2020; and Gain probability as the ratio of empty cells in 2019 that became occupied in 2020. Probabilities were averaged at 20 cm intervals along the intertidal zone occupied by mussels, ranging from 3.2 to 5.8 m above the LAT (Barbosa et al., 2022). In addition, the error of classification (1- accuracy) was used to determine the true observed dynamics and the variability of the resultant mussel distribution area, i.e., SD, at each microhabitat. Using index 1 for the Occupied state and index 0 for the Empty state, the probability of classifying state *i* as state *j* was denoted as $r_{ij}$. First, the bias of the observed loss, gain, and persistence probabilities and occupation rate due to the classification errors were demonstrated (see Supplementary material 2: "Distribution of transition probabilities: mathematical proof"). Second, the bias was corrected, and the "true" values of loss, gain, and persistence were estimated. Third, the distribution of the true probabilities were derived by randomly drawing values of Loss and Gain transition in *N* ($r_{10}$, $s_{10}$) where $s_{10}$ (resp. $s_{01}$) represents the standard deviation of $r_{10}$ (resp. $r_{01}$) (see Supplementary material 2: "Distribution of transition probabilities: mathematical proof"). A series of 1000 draws was run to compute the true proportions of empty cells, the proportion of occupied cells, the Gain probability, the Loss probability, and the Persistence probability. The obtained average true values were used for all further analyses (Supplementary material 2).

Furthermore, the true probabilities were determined (including the probability of loss, persistence, and gain) for each topographic micro-habitat (intertidal height-slope-orientation) and each intertidal height to evaluate the patterns among microhabitats and the intertidal pattern, respectively. To achieve this, cells with similar microhabitat features were aggregated (from the distribution dynamics layer), i.e., by intertidal height, slope, and orientation, and calculated loss, persistence, and
gain probabilities as indicated above. The slope was classified into five ranges, each spanning 10° from 0 to 60. Similarly, the orientation was divided from 0 to 360 every 30° (counterclockwise, with 0 representing North). A merged dataset containing mussel distribution dynamics and topographic features corresponding to each cell in the orthomosaic was created. To achieve this, the distribution dynamics layer was converted from raster to points, and this point layer was used to extract information from layers of topography (DSM), slope, and orientation (at a resolution of 20 × 20 cm). The topographic feature layers were derived from the DSM, which was downscaled from 1 to 20 cm resolution, using the raster R package (Hijmans, 2020). The 20 cm resolution of these features was selected to facilitate processing across the entire study area and avoid very localized roughness effects. Such resolution of topographic features was shown to be related to the occupation rate of mussels (Barbosa et al., 2022), supporting our decision.

To identify the main features related to true mussel distribution dynamics, we evaluated the amount of variance explained by each feature and their potential interactions by fitting linear models (LMs) in R (R Core Team, 2020). Values of true gain from the East shore and true loss, true gain, and true persistence from the West shore were transformed to normalize the distribution of the LMs residuals by applying the arcsine square root transformation. Then, the variance explained was calculated using the proportion of the total variation explained by a variable, denoted as η² Lane.

$$η^2 = SS_{\text{effect}} / SS_{\text{total}}$$

Where SS effect is the sum of squares due to an effect or factor, and

SS total is the total sum of squares. The proportion of the total variation explained by each variable (η²) was represented as a percentage in a Venn diagram for loss, persistence, and gain from each shore. Intertidal height was treated as a continuous variable represented by ranges of 20 cm along the intertidal zone occupied by mussels. The orientation topographic feature (independent variable) was treated as a categorical variable with several levels determined by a range of 30°, from 0 to 350. The slope topographic feature was treated as a continuous variable. Cells with a slope higher than 60 were excluded from the analyses due to their low representation (few cells and lack of representation of all orientation angles) in the study area. The entire process of extracting layer information and creating the complete dataset was performed in R (R Core Team, 2020) using raster package functions (Hijmans, 2020).

## 2.5. Model of the dynamics of mussel's occupied area

The model represents the dynamics of a theoretical cell defined as a square area of 1 cm², which is the same scale of observations made at a given intertidal height position on a single shore. Each cell can be in two ecological states, depending on the presence (Occupied) or absence (Empty) of mussels (Fig. 2A). With the same notations as in the previous section, qij the probability to change from state i to state j. An Occupied cell has, therefore, a probability of staying occupied for one more year (persistence probability), denoted as q11. Each Occupied cell may also become Empty with a loss probability q10. The probability of change from Empty to Occupied between two consecutive years is defined as the gain probability q01 (Fig. 2A).

<figure><img src="image_4.png" /><figcaption>Fig. 2. A) States and transition probabilities of cells for a specific intertidal height range and shore. B) Temporal dynamics of a theoretical proportion of Occupied cells.</figcaption></figure>
The model was based on several assumptions. First, occupation proportions would not change over time in the absence of an extreme event (or pulse perturbation, sensu Harris et al., 2018). Some cells would become empty, and others become occupied due to regular changes (e. g., recruitment, mortality of mussels) but the whole distribution of *Occupied/Empty* cells and the loss/gain probabilities would maintain over time and result in self-organized patterns which are stable states (Liu et al., 2014). Since there were only 2 observed spatial distributions (from the observations in 2019 and 2020), there was only one possible set of derived occupation, loss, and gain probabilities corresponding to the steady state. The differences between model predictions and observations were interpreted as the results of disturbances, which are masked in the observations. The second assumption was that the sensitivity of bathymetry levels to additional disturbances would show the same patterns as the ones that have been observed between 2019 and 2020. Therefore, resilience, persistence, and return time predicted by the model likely differ from observed values, but their patterns are consistent, which makes the comparison of predictions of shores and bathymetry levels relevant.

Based on the previous assumptions, Markov chain theory predicts that the proportions of cells in *Occupied* vs *Empty* states reach equilibrium, remaining constant over time (Wootton, 2004, Fig. 2B). Consequently, mussel distribution is expected to fluctuate annually due to mortality, while the proportions of *Occupied* vs *Empty* cells remain stable unless additional mortality occur. Properties of Markov chains allowed us to derive important features characterizing the long-term dynamics of mussel population. The questions that we address are the following: how do bathymetry levels and sites compare in terms of how long cells may remain Occupied by mussels, their Occupation rate at equilibrium, and their resilience? In summary (see Supplementary material 3: "Transactions between states using Markov Chains: Mathematical proof").

- proportion of *Occupied* cells is $\lambda = \frac{q_{01}}{q_{01}+q_{10}}$

- persistence time, defined as the maximum time during which 90% of *Occupied* cells remain occupied (90th percentile) in the absence of mortality, is equal to: $a = \frac{\log(0.1)}{\log(1-q_{10})}$

- resilience index, is a metric of the capacity to return to equilibrium after a decrease due to mortality (loss events). It is equal to $r = 1 - |1 - q_{01} - q_{10}|$. This indicator ranges between 0 and 1, facilitating comparison across systems. Values close to 1 correspond to more resilient systems, following the definition given by Viejo (2009) and Grimm and Wissel (1997) as the speed at which an ecological system returns to a reference condition after a temporary disturbance. Here, the "temporary disturbance from a reference condition" is associated with cumulative mortality over a year. To avoid confusion between disturbance from the model's equilibrium and short-term mass mortality events, generally referred to as disturbance events, we use hereafter disturbance and environmental disturbance, respectively.

- return time to equilibrium (Fig. 2B) after a disturbance is $e = \frac{\log(0.01)}{\log(1-q_{10}-q_{01})}$. This metric is related to the resilience index but measures time, as cells with a high resilience index will present a shorter return time. The return time can be site/system-specific, whereas the resilience index is a standardized metric for comparison with other systems.

Transition probabilities $q_{10}$ and $q_{01}$ likely depend on the intertidal height range and shore. There were derived from the drone image analyses by aggregating all cells from the overlapped orthomosaics at each specific intertidal height range from the West or East shore. The equations and the observed proportions of *Occupied*, *Loss*, and *Gain* cells were used to estimate the two parameters, $q_{10}$ and $q_{01}$, for each shore and intertidal height range. For each shore and intertidal height range, the average and standard deviations of Markov metrics were computed by using the distribution of transition probabilities calculated to account for the image classification errors (Table 1). Then, the properties of cell

<table><caption>**Table 1**
Overall accuracy and Kappa index indicating the performance of the supervised classification results performed for each shore and date.</caption><thead><tr><th>Rocky shore</th><th>Dataset</th><th>Overall accuracy (%)</th><th>Kappa index</th></tr></thead><tbody><tr><td rowspan="2">West</td><td>June 2019</td><td>0.98</td><td>0.93</td></tr><tr><td>June 2020</td><td>0.93</td><td>0.79</td></tr><tr><td rowspan="2">East</td><td>June 2019</td><td>0.93</td><td>0.80</td></tr><tr><td>June 2020</td><td>0.94</td><td>0.82</td></tr></tbody></table>

dynamics in these habitats were computed and compared using the same metrics. In addition, the sensitivity of Markov metrics was assessed by computing the metrics for a combination of transition probabilities spanning from the lowest to the highest values that were observed.

## 3. Results

### 3.1. Accuracy of the identified mussels' distribution dynamics

The accuracy of the image classification process was very good, with more than 79 % agreement between cells identified as mussel cells and the true mussel cells, as indicated by the Kappa index, across the two shores and dates studied (Table 1). The overall accuracy of the image classification was high, with 93 %-98 % of evaluated cells correctly classified among the processed image datasets (Table 1).

### 3.2. Mussel's distribution dynamics

All the following results represent the corrected values after accounting for the errors of transition probabilities presented above and were calculated based on the classification errors (see Methods). From the initially covered area in June 2019 (44.2 m²), 80 % (35.5 m²) was lost, and 32 % (14.3 m²) was gained on the West shore, whereas in the East (302.3 m² in 2019), 55 % (176.9 m²) of the covered area was lost, and 7 % (20.7 m²) was gained (Table 2, Fig. 3). The gain area was, thus, smaller than the loss area, resulting in a net decrease of the covered area of 20 % and 44 % on the West and East shores, respectively (Table 2, Fig. 3).

The probabilities of loss, persistence and gain varied principally between the West and East shores and among intertidal height positions, slope and orientation (Fig. 4). Between shores, the West presented higher loss probabilities and lower persistence probabilities than the East shore, mainly in the mid-intertidal area, whereas the gain probabilities were low on both shores and along the entire intertidal height gradient (Fig. 4). The loss probability and persistence probability varied from 0 to 1 in the West and East shore (Fig. 4). The gain probability varied between 0 and 0.1 in the microhabitats of the West and East shore (Fig. 4). Along the intertidal height gradient, the loss and persistence probabilities presented a pattern with lower loss and higher persistence probability in the mid-intertidal height, mainly on the East shore, but no evident pattern was observed in the gain probability on both shores (Fig. 4). Regarding the microhabitat slope and orientation, there was no evident pattern either at specific intertidal positions or shore (Fig. 4).

Analyzing the proportion of variance explained by each of these factors also evidenced the higher importance of intertidal height over the slope and orientation of microhabitats on the probabilities of loss, persistence, and gain, in both rocky shores (Fig. 5, Table S1 in Supplementary material 1). However, the topographic features effect on the two shores were different, with higher effects observed in the East compared to the West shore. The total variance not explained by either intertidal height, slope, and orientation (white area; Fig. 5) was higher in the West shore (68.2 %-69.7 %) than in the East (19.9 %-56.8 %), indicating an effect of the shore on the observed patterns of loss, persistence, and gain. A large part of the variance observed in the West shore would, therefore, be associated with another factor unrelated to the topographic features studied here. On the East shore, intertidal
<caption>**Table 2**</caption>

Summary of distribution dynamics (true values, i.e., corrected by error of image classification) in the West and East shores in ‘Le Petit Minou’. Percentages are calculated with respect to the total initial area of the corresponding shore (2019; 100 %).

<table><thead><tr><th>Shore</th><th>Occupied Initial area, 2019 (m²)</th><th>Occupied Final area, 2020 (m²)</th><th>Net area difference m² (%)</th><th>Loss surface m² (%)</th><th>Persistence surface m² (%)</th><th>Gain surface m² (%)</th></tr></thead><tbody><tr><td>West</td><td>44.2</td><td>23.0</td><td>-21.2 (47.9 %)</td><td>35.5 (80 %)</td><td>8.6 (20 %)</td><td>14.4 (32 %)</td></tr><tr><td>East</td><td>302.3</td><td>155.1</td><td>-147.2 (48.7 %)</td><td>167.9 (56 %)</td><td>134.3 (44 %)</td><td>20.7 (7 %)</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Mussel coverage distribution dynamics (±SD) between June 2019 and June 2020 in the West and East rocky shore in ‘Le Petit Minou’, France. Inset images from the West (top panels) and East shore (bottom panels) show mussel occupancy in 2019 and 2020, along with the resultant dynamics. Areas where individuals were present in 2019 but absent in 2020 (loss), areas without mussels in 2019 but with mussel presence in 2020 (gain), and areas where mussels were present in both 2019 and 2020 (persistence) are indicated. Yellow stars indicate the location of inserted images; scale bars are in meters. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** The probability of loss, persistence, and gain area for all microhabitats with different slopes (increasing from the center to the border; y axis), orientation (cardinal orientations indicated with letters; N: North, E: East, S: South, and W: West), and intertidal height (increasing from bottom to the top panels). Loss probability is the ratio of occupied cells in 2019 that become empty in 2020; Persistence probability is the ratio of occupied cells in 2019 that continue to be occupied in 2020; and Gain probability is the ratio of empty cells in 2019 that become occupied in 2020. Probabilities represent the average of true values per microhabitat after accounting for classification errors.</figcaption></figure>

height explained 75.6 % and 75.6 % of the loss and persistence variation, respectively, while it explained 23.0 % and 23.0 % on the West shore. Contrarily, the variance explained by slope and orientation was small and relatively similar between West and East shores, ranging from 0.2 to 8.3 % and from 0.2 to 2.3 % for slope and orientation, respectively.

Considering the primary effect of intertidal height on mussel distribution dynamics, we explored how mussels' loss, persistence, and gain probabilities were distributed along both East and West intertidal

heights. The distribution of their probabilities showed different patterns along the intertidal height (Fig. 6). On both shores, the persistence showed an unimodal distribution along the intertidal height, distributed on a wider range of values on the East shore and centered on 4.5 m above the LAT, relative to the West where values were observed over a lower range, between 4.3 m and 5.3 m, and centered on 4.7 m above the LAT. Gain probabilities showed a more similar unimodal pattern between shores, centered on 4.5 m above the LAT (Fig. 6; note the different scales). The loss probability presented an inverse pattern relative to the
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Venn diagrams with the percentage of explained variance, η² (based on the sum of squares, Table S1 in Supplementary material 1), that was related to each of the three topographic features evaluated (circles) and the interactions among them (overlapped areas). The percentage of variance not explained by the topographic features is indicated with the white background area. Percentage of variance in bold indicates statistically significant variables/interactions (p < 0.05).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Relative distribution dynamics are represented as the probability of loss, persistence, and gain at specific intertidal height ranges on the West and East shore. Loss probability is the ratio of occupied cells in 2019 that become empty in 2020; Persistence probability is the ratio of occupied cells in 2019 that continue to be occupied in 2020; and, Gain probability is the ratio of empty cells in 2019 that become occupied in 2020. True values of probabilities were averaged at each intertidal height range. Variability bars indicate the standard deviations (SD). The left axis represents Loss and Persistence probabilities, whereas the right axis represents Gain probability.</figcaption></figure>

persistence, with a bimodal distribution on both shores, with maxima observed at both margins through the top and bottom intertidal heights (Fig. 6).

### 3.3. Properties of the mussel dynamics predicted by markov model

Based on the primary effect of the rocky shore and intertidal height on the probabilities of loss, persistence, and gain, we estimated model predictions along the intertidal height on each rocky shore. The uncertainty in the estimation of loss and gain probabilities yielded very small standard deviations in all the metrics that we computed (Fig. 7).

Standard deviations of occupation rates were smaller than their average values.

The comparison between observation and model predictions indicated that in 2019 and 2020, the observed occupation rates were greater than the predicted occupation rates at equilibrium (Fig. 7). For instance, the maximum observed values on the East shore were around 0.20 in 2020 and 0.15 in 2019, as opposed to maximum predicted values of 0.02. This is due to the loss and gain probabilities, which were based on the changes measured between 2019 and 2020, and that we used in our model parametrization. Since there was a sharp loss probability and, thus, a decrease in the occupation rate between the 2 years, the model
<figure><img src="image_9.png" /><figcaption>Fig. 7. Occupation rate (±SD) for 2019 and 2020, as well as the predicted occupation rate at equilibrium for specific intertidal height ranges on the West and East shore.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 8. Predicted values (±SD) along the intertidal height of the West and East shores, A) Return time (years) to equilibrium; B) resilience index (no unit; values between 0 and 1); C) persistence time (years) defined as the time during which 90 % of the occupied cells remain occupied (quantile 90).</figcaption></figure>
predicts that, if this loss is maintained over the years, the Occupation proportion would remain very low. The sensitivity analysis showed that predicted metrics varied at different values of loss probabilities, which were much greater than gain probabilities (Fig. S2). In all intertidal heights and both rocky shores, the variability of the predicted metrics at different transition probabilities remained smaller than the average metric values (Fig. S2).

The time to return to stable *Occupied* and *Empty* proportions after a disturbance varied between 1 and 7 years (Fig. 8A). The resilience index was above 0.5 for all combinations of shore and intertidal height, indicating a high capacity to return to equilibrium (Fig. 8B). Resilience is inversely (but not linearly) correlated to Return time (see Methods). Resilience was higher on the West than on the East shore, which corresponded with a slightly shorter return time shown by the simulations (Fig. 8B). For both shores, resilience was lower in the mid-intertidal zone, around 4.6 m above the LAT. The persistence time and resilience index distribution presented opposite patterns along the intertidal height (Fig. 8B–C). Persistence time was higher on the East shore and varied along the intertidal height with maximum values of 3.5 years at intermediate intertidal height (Fig. 8C).

## 4. Discussion

Using drone images and a Markov chain model, we showed that the covered area by mussels (*Mytilus galloprovincialis*) on rocky shore habitats over a year was highly dynamic and spatially heterogeneous. The proportion of area with loss, gain, and persistence of mussels differed between the two rocky shores studied and among intertidal height ranges but they were quite similar among microhabitats of different slopes and orientations. Our results then support our first hypothesis about the persistence of covered areas determining the higher occupation rate (occupied proportion of the available area) at intermediate intertidal height and its decrease at the top and bottom edges. On the contrary, our second hypothesis about the effect of microhabitat slope and orientation on the dynamics of occupancy was not supported since the variation related to these features represented less than 9 % of the overall variation. The modeling approach showed a higher capacity of the population to recover from loss/mortality of mussels at the top and bottom intertidal edges (maximum resilience index greater than 0.5 and return time lower than 3 years). Furthermore, it indicated a potential of mussel areas to persist over 3.6 years at intermediate intertidal height areas in the absence of mortality events. In the following, we discuss the spatial heterogeneity of mussel dynamics and sociated predicted long-term population properties (resilience and persistence time) and how they may relate to environmental factors or population traits. Finally, we discuss the value of combining observations and modeling when studying population ecology.

### 4.1. Mussels' occupancy and temporal dynamics

The mussel-covered areas in “Le Petit Minou” decreased substantially, between 50 % and 80 %, from June 2019 to June 2020 on both rocky shores and were only partially recovered. This highly dynamic character is not specific to the studied population and has also been described for mussels in different regions. For example, monitoring of a wild population of *M. californianus* in the Northwest Pacific showed that natural mortality linked to environmental disturbance in a rocky shore led to annual losses from 30 % to 40 % during the 1970–1978 period (Paine and Levin, 1981). Other studies, based on experimental mortality (by removing individuals) approaches and monitoring of small areas on rocky shores in Rhode Island, Northwest Atlantic, highlighted highly dynamic seasonal cover loss and recovery of mussels in the middle intertidal zone (0.5–1.0 m above mean lower low water). By monitoring fixed quadrats over a 3-yr period (2001–2003), Carrington et al. (2009) showed a cover loss of about 50 % in the first year followed by stability and recovery or gain in the third year, a similar recovery time as our

prediction for the East shore from “Le Petit Minou”. Another quadrat showed about 20 % loss area in the first years, followed by a similar recovery and a subsequent loss of ca. 40 % which was not recovered in the second year, but increased in the third year reaching ca. 100 % loss, which was recovered to 70 % coverage in only ca. 4 months. These results by Carrington et al. (2009) showed the high spatial variability and temporal dynamic in loss and recovery (or gain) at a single intertidal range.

Our methodological approach, using remote-sensing imagery, provides extensive coverage of the mussel habitat and has enabled us to show that mussels' distribution dynamics varied principally along the intertidal height gradient on the two shores, with between 8 and 58 % of the variation related to it. This supports our first hypothesis that the persistence of mussel aggregations determines the higher occupancy at intermediate intertidal height. In other words, our results showed that the spatial distribution of mortality, observed as loss areas, could determine the unimodal distribution of mussels along the intertidal gradient. This suggests an optimum habitat for *M. galloprovincialis* when intertidal height is between 3.8 and 5 m above the LAT in our study area. Based on our mussels' dynamics observation, Markov chain modeling showed that the pattern of persistence time was consistent with the observed persistence probabilities and, thus, higher at this optimum intertidal height. Here mussels were predicted to reach a maximum persistence time of ca. 3.6 years at intermediate intertidal height and about 1 year at the distribution edges. That estimation appears much lower than the maximum age of mussel individuals recorded in our study area, “Le Petit Minou”, about 12 years old based on sclerochronology analyses (Barbosa et al., 2021). Such a mismatch could reflect the randomness of a few individuals reaching such an age of 12 years. Furthermore, our modeling approach simplifies the dynamics of mussels in a square area and assumes a complete change of state in such an area, whereas in reality, even in a disturbed area, some individuals can persist and create mussel aggregations of heterogeneous age as observed *in-situ* by Barbosa et al. (2021).

In addition, our modeling approach showed an intertidal pattern of the resilience index and return time, which suggested that the population distribution is highly resilient, mainly in the top and bottom edges of the intertidal range. It could sound contradictory since the areas with a lower resilience are the most occupied. However, mid-intertidal areas had a lower loss probability, which suggests that the ecological strategy is likely to increase their resilience in areas with higher probability of loss. It is important to note that a lower loss probability in the mid-intertidal zone does not necessarily indicate a smaller loss area compared to the intertidal edges. In fact, a larger loss area in the mid-intertidal, where the covered area is greater, can result in a smaller loss probability than a similar loss area at the edges where the mussel covered area is smaller, because the probability is related to a proportion. Higher resilience was also recorded at higher intertidal positions, with higher environmental stress in macroalgal communities (Viejo, 2009). The predicted return time at Le Petit Minou mid-intertidal height was in concordance with slow recovery in mussel coverage in wave-exposed areas in Patagonia, from 5 to 15 % recovery after 3 years of experimental mortality (Bertness et al., 2006). The trade-off between loss probability and resilience could explain how the species can inhabit the upper and lower intertidal areas, where, without a high recovery capacity, mussels would completely disappear over time.

Contrary to expected, the spatial distribution of loss and gain probability (i.e., the dynamics of mussels' distribution) recorded between 2019 and 2020 was only marginally related to microhabitat slope and orientation. This result may seem contradictory to observations made by Harley (2008), who found higher mortality in microhabitats oriented to the south during high-temperature periods. However, our one-year period dynamics balance must limit our ability to identify the effect of microhabitat features at specific seasons. For instance, mortality driven by environmental conditions such as elevated temperatures may predominantly affect microhabitats with specific orientations during
summer, whereas mortality resulting from wave action may be more prevalent in winter and spring, affecting microhabitats with distinct orientations. Furthermore, our results do not explain the observations made by Barbosa et al. (2022), who showed a link between the distribution of mussels on the East shore of “Le Petit Minou” in 2019 and microhabitats’ features. The occupancy pattern recorded on the East shore should, thus, result from a different process than mortality between 2019 and 2020 or even before 2019. Another potential explanation could be that a cumulative effect of longer-term processes would produce the observed occupancy pattern linked to orientation and slope, and short-term mortality would be mainly linked to the intertidal height profile and the effect of shore (i.e., mesoscale).

## 4.2. Potential factors driving the distribution dynamics

Our observations of a net decrease in the total occupation rate indicated that the gained area was not enough to compensate for the loss of individuals in a year. The system was, therefore, not at equilibrium, i. e., it does not have a constant proportion of *Occupied* vs. *Empty* area over time, which calls for an explanation based on one or more sources of mortality. The low gain area observed in the present study indicated slow recolonization during 2019–2020. Recolonization of loss areas generally occurs by the growth of individuals still in place and movement of individuals from the perimeter of the loss areas, also called clearings or gaps, while recruitment in such bare areas is almost insignificant (Paine and Levin, 1981; Littorin and Gilek, 1999). Massive recruitment is observed in other regions after mass mortality events, such as in rocky intertidal rocky-shore areas from Norway (Strohmeier et al., 2024), the Northwest Atlantic and Northeast Pacific (e.g., Scrosati, 2024; Scrosati and Cameron, 2024), and in areas where mussel seeds are collected for cultivation such as in Galicia, Spain (Pérez Camacho et al., 1991). However, such areas with evident recolonization were not observed in Le Petit Minou during the studied year. Potential recruitment may occur through the settlement of larvae in filamentous algae at lower intertidal followed by a migration of small recruits to the mid and upper intertidal, as registered in other regions (Bayne, 1964). We don’t have evidence of such a process in our study area and it was out of the scope of this study. However, evidence exists about recruitment inside the existing mussel aggregations from Le Petit Minou, which is -linearly related to mussels’ density (Barbosa et al., 2022). Another potential explanation for the observed low gain area in only one year is the slow growth rate of mussels on both shores, with less than 1 cm growth over a year (Barbosa et al., 2021). Mussel bed mortality is mainly related to environmental disturbance events of wave action and temperature (Harger and Landenberger, 1971; Paine and Levin, 1981; Denny, 1987; Harley and Helmuth, 2003; Guichard et al., 2003; Harley, 2008; Carrington et al., 2009; Seuront et al., 2019). Combined conditions of high-temperature stress and wave action cause mass mortality by dislodgement in summer, when mussels’ attachment is generally weaker (Carrington, 2002; Carrington et al., 2008, 2009). The high loss observed in our study area could result from this type of extreme event. Indeed, in our study region, an event combining high spring and summer temperatures compared to long-term records and a period of high waves occurred in early August (Fig. S3 in Supplementary material 1). The occurrence and intensity of marine heatwaves have been increasing over the last decade because of global warming in the Bay of Biscay, English Channel, and Bay of Brest (Simon et al., 2023). Consequently, environmental disturbances of mussel distribution related to marine heatwaves are likely to occur at a larger spatial scale and a higher frequency.

Sand dynamics and predation can also affect mussel distribution, mainly at the lower intertidal level (Peterson, 1979; Rilov and Schiel, 2006; Zardi et al., 2008; Bagur et al., 2022). Sand dynamics increases mussels’ mortality and determines mussels’ lower limit of distribution in other areas (Zardi et al., 2008; Bagur et al., 2022). Sand intertidal level in the studied site was observed to change by about 1 m in height during the winter of 2019 (Feurer et al., 2024). These movements of the sand

and changes to the site’s overall topography may have caused additional environmental disturbance, helping to explain the low rate of gain observed. Sea star predators are not the main predators in our study area (Tsuchiya and Retiere, 1988). Other predators such as whelks (Nucella lapillus), crabs (Carcinus maena), and seagulls can affect the mussel distribution in the Northeast Atlantic (Dudgeon and Petraitis, 2019; Rilov and Schiel, 2006), but there is no quantitative information available about the predation in our study area. Still, whelks may be one of the main predators due to their high abundance and main distribution along the mid- and low-intertidal in the area (Fred Jean personal communication). On the other hand, there is a decrease in crabs’ abundance on wave-exposed rocky shores (Dudgeon and Petraitis, 2019; Silva et al., 2020), which mainly inhabit lower intertidal or subtidal, and go up and down depending on the tide (Naylor, 1962; Silva et al., 2020). The crabs may regulate mussels or whelk abundance, although crab abundance may be medium to low in our site due to its high wave exposure. This would not play a main role in the limit of mussel range, as also observed in Patagonia, where sand dynamics limit mussels’ lower distribution (Bertness et al., 2006). We speculate that predation might contribute to mid and low intertidal mortality and potentially prevent recruitment in bare areas (Rilov and Schiel, 2006), while sand dynamics could set the lowest distribution limit.

The spatial variability in mortality sources and recruitment combined with differences in the biological features of mussel aggregations between shores, could explain the observed differences in mussel cover dynamics between the East and West shores. Based on our results, mussel aggregations in the West appear less stable, presenting a higher probability of loss and consequently lower persistence than those in the East. Due to the distinct orientation of both shores, differences in wave exposure and sun exposure are expected between the East and West shores, with consequently differences in the average temperatures experienced by individual mussels (Helmuth and Denny, 2003; Gilman et al., 2006; Seabra et al., 2011). Potentially higher sun exposure could intensify the mussels’ temperature and, thus, the combined effect of waves and thermal stress. Furthermore, patterns of mussel occupancy were related to mesoscale orientation differences between shores (Barbosa et al., 2022). Population features, such as density and crowding structure, could also influence the effect of environmental disturbance events. Environmental disturbance would spread through direct neighbors, and the formation of hummocks tends to decrease the attached surface when density increases (Tanaka and Magalhães, 2002; Guichard et al., 2003; Gutiérrez et al., 2015). In reality, the spreading of the environmental disturbance among neighbors causes a more significant loss than the one produced by the actual disturbance (Denny, 1987; Guichard et al., 2003), giving density-linked stochasticity an important role in the long-term dynamic of mussels (Wootton and Forester, 2013). That is, the stochasticity of disturbance events can be combined with density-linked processes such as dislodgement. Then, the extent of loss area would result from the combination of stochasticity due to environmental events at the microhabitat scale and the dislodgement spread dependent on the spatial heterogeneity of mussels’ density. Similarly, the density of aggregations and substrate type can influence the recolonization, i.e., gain, of mussels (Bertolini et al., 2018). The density of mussels is higher in the West than in the East shore, with a higher recruitment and lower growth rate (Barbosa et al., 2021). All these features could explain the higher resilience estimated for the West compared to the East shore.

The major heterogeneity in mussel dynamics along the intertidal height, compared to the microhabitat slope and orientation, could reflect the marked tidal dynamics throughout the year compared to changes in horizontal spatial conditions. A tidal cycle would be relatively similar throughout the entire year (e.g., the lower intertidal remains always submerged for a longer period), whereas no particular pattern favoring specific orientation and slope emerges due to the high seasonality of temperature and waves exposure at a small scale (Littorin and Gilek, 1999; Gilek et al., 2001). For instance, Helmuth and Denny
(2003) highlighted the difficulty of classifying microhabitats based on the hydrodynamic forces experienced by individuals living there due to the high temporal and spatial heterogeneity in wave microscale forces (see also Burel et al., 2019; Focht and Shima, 2020; Gaylord, 1999).

Other processes, such as the spreading of environmental disturbance and density-linked stochasticity, as well as some limitations of our methods, could create noise that impacts the identification of a micro-habitat effect on the probability of loss. Density-linked stochasticity discussed above as a potential cause of differences between shores would also play a role in the differences among microhabitats. Then, we do not discard the potential relationship between slope and orientation on the probability of an initial loss of individuals, which could then spread to other microhabitats. Between the limitations of our methods, the resolution of microhabitat features (slope and orientation; 20 × 20 cm resolution), could affect the analyses.

### 4.3. Interest and limits of the modeling approach

Markov chain theory allowed us to derive resilience, return time to equilibrium, persistence time, and occupation rate at equilibrium from the transition matrix between two states, e.g., Empty cells or cells Occupied by mussels. Though we only gave the equations for this simple system, the same mathematical steps would be easily extended for systems with more than 2 states. The approach relies upon the estimation of transition probabilities, that is, the probability of loss and gain for the sites we monitored. Parameter estimation is usually based on monitoring of quadrats (Tsujino et al., 2010) or experiments (Guichard et al., 2003; Wootton, 2004) and, the accuracy of estimated transition probabilities improves with the length of the time series. In our study, we made observations on only two different dates, and our survey showed an important change in the areas covered by mussels. This change reflects the existence of a disturbance that decreased the occupation proportions from 2019 to 2020. We chose to use the observed loss and gain probabilities to approximate the characteristics of the long-term dynamics of mussel occupation. Therefore, our predictions of the occupation rate at equilibrium likely underestimate the real values that would be more accurately assessed with long-term monitoring (e.g., Pandolfo et al., 2017). However, the sensitivity analysis showed that the variability of the metrics due to the uncertainty of transition probabilities remained much smaller than the average predicted metric values. Observations over a longer period of time would presumably yield the same patterns linked to the effects of intertidal height and shore orientations, even though the values of predicted persistence time, resilience, and return time would differ from our actual results. The patterns observed could only change if the probability of loss, gain, or persistence at each intertidal height varies at different proportions over the years. This would mean, for instance, that predation pressure should present a different spatial pattern depending on the year, which is less probable.

Another source of error was linked to the uncertainty of occupation rate and transition probabilities based on the complete coverage of the area on one hand and the estimation of classification errors on the other hand. Our calculations prove that the uncertainty of the parameters of interest was negligible but allowed us to correct for a slight bias that was further considered in our Markov model. This mathematically based errors of transitions can improve the studies of temporal dynamics based on aerial images classifications, which are generally ignored when analyzing populations resilience and persistence (Olofsson et al., 2014).

## 5. Conclusions and perspectives

Our analyses allowed for an evaluation of mussel distribution dynamics over one-year, with changes observed at a spatial resolution of 1 cm along two entire rocky shores. Our results outline the role of intertidal height on the temporal and spatial dynamics of mussel distribution, the occurrence of major mortality, and the probability of long-term persistence. This work highlights the power of image analyses in the

study of spatial ecological processes related to mortality by environmental disturbances, which has been predominantly explored through experiments conducted in limited, specific areas (generally small square areas). A monitoring strategy incorporating more frequent surveys alongside ongoing temperature monitoring projects (e.g., Seabra et al., 2015; Seuront et al., 2019), would facilitate the assessment of seasonal patterns in population dynamics, capture responses to events like storms, and help identify the key drivers involved. Using drone images in conjunction with modeling approaches would be particularly apt for long-term monitoring, shedding light on the long-term persistence of populations facing changes in the frequency and intensity of marine heatwaves (Simon et al., 2023).

## CRediT authorship contribution statement

**Romina Vanessa Barbosa:** Writing – review & editing, Writing – original draft, Visualization, Methodology, Formal analysis, Data curation, Conceptualization. **Cedric Bacher:** Writing – review & editing, Visualization, Methodology, Funding acquisition, Formal analysis, Conceptualization. **Fred Jean:** Writing – review & editing, Supervision, Funding acquisition, Conceptualization. **Marion Jaud:** Writing – review & editing, Supervision, Software, Methodology, Data curation. **Jérôme Ammann:** Writing – review & editing, Software, Methodology, Data curation. **Yoann Thomas:** Writing – review & editing, Supervision, Project administration, Methodology, Funding acquisition, Conceptualization.

## Funding

This study and RVB benefited from the financial support of the Region Bretagne and Université de Bretagne Occidentale (UBO), France. The scientific committee of IUEM (France) funded student's training related to the project.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.marenvres.2025.107321.

## Data availability

Data will be made available on request.

## References

* Agisoft, 2018. Agisoft PhotoScan (Version 1.4.0) [Computer software]. Agisoft LLC. https://www.agisoft.com/.

* Arribas, L.P., Donnarumma, L., Palomo, M.G., Scrosati, R.A., 2014. Intertidal mussels as ecosystem engineers: their associated invertebrate biodiversity under contrasting wave exposures. Mar. Biodivers. 44, 203–211.

* Bagur, M., Gutiérrez, J.L., González, J.A., Arribas, L.P., Palomo, M.G., 2022. Physical rather than biotic factors set the lower limit of mussel beds in a horizontal rocky intertidal platform. J. Exp. Mar. Biol. Ecol. 548, 151680.

* Barbosa, R.V., Bacher, C., Jean, F., Thomas, Y., 2021. Individual Traits and Population Parameters of a rocky-shore Mussel Population.

* Barbosa, R.V., Jaud, M., Bacher, C., Kerjean, Y., Jean, F., Ammann, J., Thomas, Y., 2022. High-resolution drone images show that the distribution of mussels depends on microhabitat features of intertidal rocky shores. Remote Sens. 14, 5441.

* Bayne, 1964. Primary and secondary settlement in Mytilus edulis. L. (Mollusca). J. Anim. Ecol 513–523.

* Bertness, M.D., Crain, C.M., Silliman, B.R., Bazterrica, M.C., Reyna, M.V., Hildago, F., Farina, J.K., 2006. The community structure of Western Atlantic patagonian rocky shores. Ecol. Monogr. 76, 439–460.
Bierne, N., Borsa, P., Daguin, C., Jollivet, D., Viard, F., Bonhomme, F., David, P., 2003. Introgression patterns in the mosaic hybrid zone between *Mytilus edulis* and *M. galloprovincialis*. Mol. Ecol. 12, 447–461.

Bertolini, C., Montgomery, W.I., O'Connor, N.E., 2018. Habitat with small inter-structural spaces promotes mussel survival and reef generation. Mar Biol 165, 163.

Borthagaray, A.I., Carranza, A., 2007. Mussels as ecosystem engineers: their contribution to species richness in a rocky Littoral community. Acta Oecol. 31, 243–250.

Burel, T., Schaal, G., Grall, J., Le Duff, M., Chapalain, G., Schmitt, B., Gemin, M., Boucher, O., Ar Gall, E., 2019. Small-scale effects of hydrodynamics on the structure of intertidal macroalgal communities: a novel approach. Estuar. Coast Shelf Sci. 226, 106290.

Carrington, E., 2002. Seasonal variation in the attachment strength of blue mussels: causes and consequences. Limnol. Oceanogr. 47, 1723–1733.

Carrington, E., Moeser, G.M., Dimond, J., Mello, J.J., Boller, M.L., 2009. Seasonal disturbance to mussel beds: field test of a mechanistic model predicting wave dislodgment. Limnol. Oceanogr. 54, 978–986.

Carrington, E., Moeser, G.M., Thompson, S.B., Coutts, L.C., Craig, C.A., 2008. Mussel attachment on rocky shores: the effect of flow on byssus production. Integr. Comp. Biol. 48, 801–807.

Choi, F., Gouhier, T., Lima, F., Rilov, G., Seabra, R., Helmuth, B., 2019. Mapping physiology: biophysical mechanisms define scales of climate change impacts. Conserv Physiol 7, coz028.

Cohen, J., 1960. A coefficient of agreement for nominal scales. Educ. Psychol. Meas. 20, 37–46.

Congalton, R.G., 1991. A review of assessing the accuracy of classifications of remotely sensed data. Remote Sens. Environ. 37, 35–46.

Curd, A., Chevalier, M., Vasquez, M., Boyé, A., Firth, L.B., Marzloff, M.P., Bricheno, L.M., Burrows, M.T., Bush, L.E., Cordier, C., Davies, A.J., Green, J.A.M., Hawkins, S.J., Lima, F.P., Meneghesso, C., Mieszowska, N., Seabra, R., Dubois, S.F., 2023.

Applying landscape metrics to species distribution model predictions to characterize internal range structure and associated changes. Glob. Change Biol. 29, 631–647.

Dayton, P.K., 1971. Competition, disturbance, and community organization: the provision and subsequent utilization of space in a rocky intertidal community. Ecol. Monogr. 41, 351–389.

Denny, M.W., 1987. Lift as a mechanism of patch initiation in mussel beds. J. Exp. Mar. Biol. Ecol. 113, 231–245.

Dudgeon, S.R., Petraitis, P.S., 2019. Rocky intertidal shores of the north-west Atlantic Ocean. Chapter. In: Hawkins, Stephen J., Bohn, Katrin, Firth, Louise B., Williams, Gray A. (Eds.), Interactions in the Marine Benthos: Global Patterns and Processes, Systematics Association Special Volume Series. Cambridge University Press, Cambridge, pp. 61–89.

Feurer, D., Bemis, S., Couloma, G., Mabrouk, H., Massuel, S., Barbosa, R.V., Thomas, Y., Ammann, J., Vinatier, F., 2024. Broader applications of the Time-SIFT method. In: Mallet, C. (Ed.), Multitemporal Earth Observations Image Analysis. https://doi.org/10.10002/9781394306657.ch1.

Feurer, D., Vinatier, F., 2018. Joining multi-epoch archival aerial images in a single SfM block allows 3-D change detection with almost exclusively image information. ISPRG J Photogrammetry Remote Sens 146, 495–506.

Focht, R.C., Shima, J.S., 2020. Acceleration loggers reveal fine-scale heterogeneity in wave exposure along an open Coast. Estuar. Coast Shelf Sci. 233, 106507.

Garza, C., 2019. Landscape ecology in the rocky intertidal: opportunities for advancing discovery and innovation in intertidal research. Curr Landscape Ecol Rep 4, 83–90.

Gaylord, B., 1999. Detailing agents of physical disturbance: wave-induced velocities and accelerations on a rocky shore. J. Exp. Mar. Biol. Ecol. 239, 85–124.

Gilek, M., Littorin, B., Saetre, P., 2001. Spatial patterns of abundance and growth of *Mytilus edulis* on boulders in the northern Baltic sea proper. Hydrobiologia 452, 59–68.

Gilman, S.E., Wethey, D.S., Helmuth, B., 2006. Variation in the sensitivity of organismal body temperature to climate change over local and geographic scales. Proc. Natl. Acad. Sci. 103, 9560–9565.

Gomes, I., Teixeiro, L., Bueno-Pardo, J., Albuquerque, R., Pérez-Jorge, S., Oliveira, E.R., Alves, F.L., Queiroga, H., 2018. What's a picture really worth? On the use of drone aerial imagery to estimate intertidal rocky shore mussel demographic parameters. Estuar. Coast Shelf Sci. 213, 185–198.

Grimm, V., Wissel, C., 1997. Babel, or the ecological stability discussions: an inventory and analysis of terminology and a guide for avoiding confusion. Oecologia 109, 323–334.

Guichard, F., Bourget, E., Agnard, J.-P., 2000. High-resolution remote sensing of intertidal ecosystems: a low-cost technique to link scale-dependent patterns and processes. Limnol. Oceanogr. 45, 328–338.

Guichard, F., Bourget, E., Robert, J., 2001. Scaling the influence of topographic heterogeneity on intertidal benthic communities: alternate trajectories mediated by hydrodynamics and shading. Mar. Ecol. Prog. Ser. 217, 27–41.

Guichard, F., Halpin, P.M., Allison, G.W., Lübchenco, J., Menge, B.A., 2003. Mussel disturbance dynamics: signatures of oceanographic forcing from local interactions. Am. Nat. 161, 889–904.

Gutiérrez, J., Palomo, M., Bagur, M., Arribas, L., Soria, S., 2015. Wave action limits crowding in an intertidal mussel. Mar. Ecol. Prog. Ser. 518, 153–163.

Harger, J.R.E., Landenberger, D.E., 1971. The effect of storms as a density dependent mortality factor on population of sea mussels. Veliger 14, 195–201.

Harley, C., 2008. Tidal dynamics, topographic orientation, and temperature-mediated mass mortalities on rocky shores. Mar. Ecol. Prog. Ser. 371, 37–46.

Harley, C.D.G., Helmuth, B.S.T., 2003. Local- and regional-scale effects of wave exposure, thermal stress, and absolute versus effective shore level on patterns of intertidal zonation. Limnol. Oceanogr. 48, 1498–1508.

Harris, R.M.B., Beaumont, L.J., Vance, T.R., Tozer, C.R., Remenyi, T.A., Perkins-Kirkpatrick, S.E., Mitchell, P.J., Nicotra, A.B., McGregor, S., Andrew, N.R., Letnic, M., Kearney, M.R., Wernberg, T., Hutley, L.B., Chambers, L.E., Fletcher, M.-S., Keatley, M.R., Woodward, C.A., Williamson, G., Duke, N.C., Bowman, D.M.J.S., 2018. Biological responses to the press and pulse of climate trends and extreme events. Nature Clim Change 8, 579–587. https://doi.org/10.1038/s41558-018-0187-9.

Helmuth, B., 2002. How do we measure the environment? Linking intertidal thermal physiology and ecology through Biophysics1. Integr. Comp. Biol. 42, 837–845.

Helmuth, B., Choi, F., Matzelle, A., Torossian, J.L., Morello, S.L., Mislan, K.A.s., Yamane, L., Strickland, D., Szathmary, P.L., Gilman, S.E., Tockstein, A., Hilbish, T.J., Burrows, M.T., Power, A.M., Gosling, E., Mieszowska, N., Harley, C.D.G., Nishizaki, M., Carrington, E., Menge, B., Petes, L., Foley, M.M., Johnson, A., Poole, M., Noble, M.M., Richmond, E.L., Robart, M., Robinson, J., Sapp, J., Sones, J., Broitman, B.R., Denny, M.W., Mach, K.J., Miller, L.P., O'Donnell, M., Ross, P., Hofmann, G.E., Zippay, M., Blanchette, C., Macfarlan, J.A., Carpizo-Ituarte, E., Ruttenberg, B., Peña Mejía, C.E., McQuaid, C.D., Lathlean, J., Monaco, C.J., Nicastro, K.R., Zardi, G., 2016. Long-term, high frequency in situ measurements of intertidal mussel bed temperatures using biomimetic sensors. Sci. Data 3, 160087.

Helmuth, B., Denny, M.W., 2003. Predicting wave exposure in the rocky intertidal zone: do bigger waves always lead to larger forces? Limnol. Oceanogr. 48, 1338–1345.

Helmuth, B., Yamane, L., Lalwani, S., Matzelle, A., Tockstein, A., Gao, N., 2011. Hidden signals of climate change in intertidal ecosystems: what (not) to expect when you are expecting. J. Exp. Mar. Biol. Ecol. 400, 191–199.

Hijmans, D.M., 2020. raster: Geographic Data Analysis and Modeling. R package version 3.6-32. https://github.com/rspatial/raster.

Hill, M.F., Witman, J.D., Caswell, H., 2002. Spatiotemporal variation in markov chain models of subtidal community succession. Ecol. Lett. 5, 665–675.

Hunt, H., Scheibling, R., 2001. Predicting wave dislodgment of mussels: variation in attachment strength with body size, habitat, and season. Mar. Ecol. Prog. Ser. 213, 157–164.

Kruskal, W.H., Wallis, W.A., 1952. Use of Ranks in One-Criterion Variance Analysis, vol. 47, pp. 583–621.

Lane, D.M., n.d. Online Statistics Education: a Multimedia Course of Study. http://onlinestatbook.com/ Project Leader: David M. Lane, Rice University.

Littorin, B., Gilek, M., 1999. A photographic study of the recolonization of cleared patches in a dense population of *Mytilus edulis* in the northern Baltic proper. Hydrobiologia 393, 211–219.

Liu, Q., Herman, M.P.J., Mooij, W.M., Huisman, J., Scheffer, H., Olff, H., van de Koppel, J., 2014. Pattern formation at multiple spatial scales drives the resilience of mussel bed ecosystems. Nat. Commun. 5 (5234). https://doi.org/10.1038/ncomms6234.

Meager, J., Schlacher, T., Green, M., 2011. Topographic complexity and landscape -temperature patterns create a dynamic habitat structure on a rocky intertidal shore. Mar. Ecol. Prog. Ser. 428, 1–12.

Naylor, E., 1962. Seasonal changes in a population of Carcinus maenas (L.) in the Littoral Zone. J. Anim. Ecol. 31 (3), 601–609. https://doi.org/10.2307/2055.

Olofsson, P., Foody, G.M., Herold, M., Stehman, S.V., Woodcock, C.E., Wulder, M.A., 2014. Good practices for estimating area and assessing accuracy of land change. Rem. Sens. Environ. 148, 42–57. https://doi.org/10.1016/j.rse.2014.02.015.

Paine, R.T., 1974. Intertidal community structure. Oecologia 15, 93–120.

Paine, R.T., Levin, S.A., 1981. Intertidal landscapes: disturbance and the dynamics of pattern. Ecol. Monogr. 51, 145–178.

Pandolfi, T.J., Kwak, T.J., Cope, W.G., Heise, R.J., Nichols, R.B., Pacifici, K., 2017. Declining occurrence and low colonization probability in freshwater Mollusk Biology and Conservation 20, 13–19.

Pérez Camacho, A., González, R., Fuentes, J., 1991. Mussel culture in Galicia (N.W. Spain). Aquaculture 94, 263–278. https://doi.org/10.1016/0044-8486(91)90122-N.

Peterson, C.H., 1979. The importance of predation and competition in organizing the intertidal epifaunal communities of Barnegat Inlet, New Jersey. Aecologia 39, 1–24. https://doi.org/10.1007/BF00345993.

Polozanska, E.S., Brown, C.J., Sydeman, W.J., Kiessling, W., Schoeman, D.S., Moore, P. J., Brander, K., Bruno, J.F., Buckley, L.B., Burrows, M.T., Duarte, C.M., Halpern, B.S., Holding, J., Kappel, C.V., O'Connor, M.I., Pandolfi, J.M., Parmesan, C., Schwing, F., Thompson, S.A., Richardson, A.J., 2013. Global imprint of climate change on marine life. Nat. Clim. Change 3, 919–925.

QGIS Development Team, 2020. QGIS Geographic Information System. Open Source Geospatial Foundation. URL http://qgis.org.

R Core Team, 2020. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/.

Rilov, G., Schiel, D.R., 2006. Thropic linkages across seascapes: subtidal predators limit effecti mussel recruitment in rocky intertidal communities. Mar. Ecol. Progr. Ser. 327, 83–93. https://doi.org/10.3354/meps327083.

Scrosati, R.A., 2024. Recovery of intertidal mussel stands three years after the severe 2021 heatwave in British Columbia, Canada. Diversity 16, 396. https://doi.org/10.3390/d16070396.

Scrosati, R.A., Arribas, L.P., Donnarumma, L., 2021. The hidden biodiversity of intertidal mussel beds under contrasting wave exposures. Bull. the Ecological Soc. Am. 102, e01811.

Scrosati, R.A., Cameron, N.M., 2024. Recolonization of intertidal mussels in Nova Scotia (canada) after their mass disappearance following the severe 2023 winter cold snap. Diversity 16, 503. https://doi.org/10.3390/d16080503.

Seabra, R., Wethey, D.S., Santos, A.M., Lima, F.P., 2011. Side matters: microhabitat influence on intertidal heat stress over a large geographical scale. J. Exp. Mar. Biol. Ecol. 400, 200–208.
Seabra, R., Whetthey, D., Santos, A., Lima, F., 2015. Understanding complex biogeographic responses to climate change. Sci. Rep. 5, 12930. https://doi.org/10.1038/srep12930.

Seuront, L., Nicastro, K.R., Zardi, G.I., Goberville, E., 2019. Decreased thermal tolerance under recurrent heat stress conditions explains summer mass mortality of the blue mussel *Mytilus edulis*. Sci. Rep. 9, 17498.

Silliman, B.R., Bertness, M.D., Altieri, A.H., Griffin, J.N., Bazterrica, M.C., Hidalgo, F.J., Crain, C.M., Reyna, M.V., 2011. Whole-community facilitation regulates biodiversity on patagonian rocky shores. PLoS One 6, e24502.

Silva, A.C.F., Hawkins, S.J., Boaventura, D.M., Brewster, E., 2020. Use of the Intertidal Zone by Mobile Predators: Influence of Wave Exposure, Tidal Phase and Elevation on Abundance and Diet. Mar. Ecol. Progr. Ser. 406, 197–210. http://www.jstor.org/stable/24873912.

Simon, A., Arbiol, C., Nielsen, E.E., Couteau, J., Sussarellu, R., Burgeot, T., Bernard, I., Coolen, J.W.P., Lamy, J.-B., Robert, S., Skazina, M., Strelkov, P., Queiroga, H., Cancio, I., Welch, J.J., Viard, F., Bierne, N., 2020. Replicated anthropogenic hybridisations reveal parallel patterns of admixture in marine mussels. Evolutionary Applications 13, 575–599.

Simon, A., Poppeschi, C., Plecha, S., Charria, G., Russo, A., 2023. Coastal and regional marine heatwaves and cold spells in the northeastern Atlantic. Ocean Sci. 19, 1339–1355.

Strohmeier, T., Strand, Ø., Haugland, B.T., Agüera, A., 2024. Baseline for the northeast Atlantic (58–70°N) intertidal mytilus species complex (mytilus spp.). Ecol. Evol. 14, e70197. https://doi.org/10.1002/ece3.70197.

Tanaka, M.O., Magalhães, C.A., 2002. Edge effects and succession dynamics in brachidontes mussel beds. Mar. Ecol. Prog. Ser. 237, 151–158.

Thomas, Y., Razafimahefa, N.R., Ménesguen, A., Bacher, C., 2020. Multi-scale interaction processes modulate the population response of a benthic species to global warming. Ecol. Model. 436, 109295.

Tsujino, M., Hori, M., Okuda, T., Nakaoka, M., Yamamoto, T., Noda, T., 2010. Distance decay of community dynamics in rocky intertidal sessile assemblages evaluated by transition matrix models. Popul. Ecol. 52, 171–180.

Tsuchiya, M., Retiere, C., 1988. Zonation of intertidal organisms and community structure of small animals associated with patches of the mussel *Mytilus edulis* L. along the rocky Coast of dinard, brittany, France. Bull. Coll. Sci. Univ. Ryukyus 54, 47–81.

Viejo, R.M., 2009. Resilience in intertidal rocky shore assemblages across the stress gradient created by emersion times. Mar. Ecol. Prog. Ser. 390, 55–65.

Wootton, J.T., 2001. Local interactions predict large-scale pattern in empirically derived cellular automata. Nature 413, 841–844. https://doi.org/10.1038/35101595.

Wootton, J., 2004. Markov chain models predict the consequences of experimental extinctions. Ecol. Lett. 7, 653–660.

Wootton, J.T., Forester, J.D., 2013. Complex population dynamics in mussels arising from density-linked stochasticity. PLoS One 8, e75700.

Zardi, G.I., Nicastro, K.R., McQuaid, C.D., Erlandsson, J., 2008. Sand and wave induced mortality in invasive (Mytilus galloprovincialis) and Indigenous (perna Perna) mussels. Mar Biol 153, 853–858.