Characterisation of extreme events waves in marine ecosystems: the case of Mediterranean Sea

We propose a new method to identify and characterise the occurrence of prolonged extreme events in marine ecosystems on the basin scale. There is a growing interest about events that can affect ecosystem functions and services in a changing climate. Our method identifies extreme events as peak occurrences over 99th percentile thresholds computed from local time series and defines an Extreme Events Wave (EEW) as a connected region including these events. The EEWs are characterised by a set of novel indexes, referred to initiation, extent, duration and strength. The indexes, associated to the 10 areas covered by each EEW, are then statistically analysed to highlight the main features of the EEWs on the considered domain. We applied the method to the winter-spring daily chlorophyll field of a validated multidecadal hindcast provided by a coupled hydrodynamic-biogeochemical model of the Mediterranean open-sea ecosystem, with 1/12° horizontal resolution. This allowed to identify the maxima of chlorophyll as exceptionally high and prolonged “blooms” and to characterise their phenomenology in the period 1994-2012. A fuzzy k-means cluster analysis on the EEWs indexes provided a bio15 regionalisation of the Mediterranean Sea associated to the occurrence of chlorophyll EEWs with different regimes.


Introduction
Extreme events affecting the Earth System have been widely investigated in hydrology and atmospheric sciences for several decades (e.g., Delaunay, 1988;Katz, 1999;Luterbacher et al., 2004;Allan and Soden, 2008;Perkins and Alexander, 2013;Tramblay et al., 2014). The study of extremes in the ocean has been concentrated mainly on the sea level (e.g., Zhang and 20 Sheng, 2015), especially in relationship to hydrology (e.g., Walsh et al., 2012), but recently there has been interest also in extreme wave height (e.g., Hansom et al., 2014), current velocity (e.g., Green and Stigebrandt, 2003) and marine heat waves (e.g., Hobday et al., 2016). However, extreme events in biosphere properties like the marine biogeochemical concentrations have received relatively little attention in the past years, despite the related heavy impacts on marine ecosystems functions and services (e.g., Zhang et al., 2010), with cascading effects at larger scales on the biogeochemical cycles (e.g., Doney, 25 2010).
Global warming, increase of atmospheric CO2 and anthropogenic eutrophication are among the major stressors of the marine ecosystems (Hoegh-Guldberg et al., 2018) and potential drivers for extreme events. Variables like SST, seawater pH, https://doi.org/10.5194/bg-2020-34 Preprint. Discussion started: 27 February 2020 c Author(s) 2020. CC BY 4.0 License. dissolved oxygen and saturation state of CaCO3 minerals are thus used as probes to monitor the marine ecosystems health (e.g., Belkin, 2009;Andersson et al., 2011;Paulmier et al., 2011). 30 In particular, some studies on marine heat waves (Hobday et al., 2016), hypoxia events (e.g. Conley et al., 2009) and low saturation state of aragonite (Hauri et al., 2013) identify extreme events from the values of a specific ecosystem variable in a site which are above/below a certain threshold, for a finite time duration. Despite there is not a common definition of "extreme event" in this context, its main emerging features result in: (i) the intensity (i.e., the absolute difference of the variable value with respect to the threshold), as a large deviation from a reference ecological state; (ii) the duration, as a 35 further stress factor on the ecosystem, eventually combined along with the intensity in an overall "severity" index (as in Hauri et al., 2013); (iii) the local character, linked to the heterogeneity of the ecosystem within the area of interest and/or due to a sparse data sampling (e.g. fixed stations). In fact, the spatial extension of the extreme event is possibly evaluated afterwards (e.g., Rabalais et al., 2002;Galli et al., 2017).
We looked for a general method able to capture the features of an "extreme event" as listed before, but also to account for 40 the persistence of the event on a certain impacted area and time duration (as Andreadis et al., 2005 andSheffield et al., 2009 proceeded in case of droughts), up to the basin scale.
The use of numerical models proves to be necessary to conduct such a study, which requires seamless and long sampling times at high frequency on the basin scale. In fact, remote sensing observations are limited by the cloud coverage and L4 data are based on filtered reconstructions (using climatology or EOF method, as in Volpe et al., 2018), which can partly hide 45 the occurrence of extremes. On the other hand, in situ-measurements do not provide the suitable spatial and temporal sampling at the basin scale and can lack in standardisation. On the contrary, numerical models provide data with continuity at high frequency in time and space. Moreover, models can account for physical and biological processes occurring in the marine ecosystems also in the sub-surface layers (e.g. vertical mixing, nutrients transport), allowing a more complete reconstruction of the dynamics of the extremes. 50 For our investigation we used the dataset provided by the MITgcm-BFM hydrodynamic-biogeochemical model of the Mediterranean Sea ecosystem at 1/12° of horizontal resolution (Di Biagio et al., 2019).
In particular, we applied the proposed method to the surface chlorophyll concentration, as Essential Ocean Variable (EOV, e.g., Muller-Karger et al. 2018), representative of the marine ecosystem state and evolution. This allowed to identify maxima of phytoplankton blooms (Desmit et al., 2018), but also positive anomalies with values too low to be actually considered 55 "blooms". In fact, due to the general oligotrophy of the basin, marked increases of phytoplankton chlorophyll are strictly considered "blooms" only in some regions (north-western Mediterranean Sea, Alboran Sea, Catalan-northern Balearic area, isolated coastal areas near some river mouths, see Siokou-Frangou et al., 2010). Our method is instead formulated to identify extreme values of chlorophyll as peaks over a threshold defined in the time series of all the basin points, i.e. at local and statistical perspective. 60 https://doi.org/10.5194/bg-2020-34 Preprint. Discussion started: 27 February 2020 c Author(s) 2020. CC BY 4.0 License.
We focused our investigation to the open-sea domain, i.e. we neglected the areas which are more directly affected by bottom and riverine dynamics (e.g., Oubelkheir et al., 2014).
The article is structured as follows. The material and methods section (Sect. 2), is divided in two parts, presenting the method to identify the extremes and the used model-derived dataset (Sect. includes a discussion about the method and the main results and it is followed by the conclusions in Sect. 5.

The method for spatio-temporal extremes investigation
The method here illustrated allows to identify, characterise and classify the extreme events in marine ecosystems, accounting for time duration and spatial extension of the events, for an ecosystem variable expressed as a concentration C(x,t). 75 Hereafter, we refer for simplicity to a daily sampling of the variable time series and to a two-dimensional variable C(x,y,t), which may be a surface variable or a quantity averaged or integrated in the vertical direction. However, the method can be easily extended to any regular time discretization and to the 3D spatial case, with few modifications in the definition of the indexes, as discussed in Sect. 4. 80

Identification
We define "extreme events" (for brevity: "extremes") the occurrences of values C(x,y,t) that are higher than a reference percentile threshold (e.g., Asch et al., 2019), computed on the whole time series of the variable. In particular, we search for the Peaks Over Threshold (POTs) referred to the 99th percentile of the time series.
Then, we identify an Extreme Events Wave (EEW) as a set of extremes that are connected in space and time. Thus, an EEW 85 tracks anomalous events that are not merely local, but that co-occur in more than a grid point and possibly are transported in space and evolve in time.
Operationally, the spatial and temporal occurrence of all the POTs can be mapped on a binary 3D matrix, representing the (2D map x day) flags of the extremes, equal to 1 for the (x,y,t) points of POTs occurrence and 0 for the points without POTs occurrence. EEW are then defined as sets of POTs occurrences that are "the closest neighbours" in space and in time. 90 The EEW definition is thus a filter on the spatio-temporal dataset: it allows to identify single events that affect a portion of the domain for a certain time period, in which all the involved points display extreme values of the variable selected. In the https://doi.org/10.5194/bg-2020-34 Preprint. Discussion started: 27 February 2020 c Author(s) 2020. CC BY 4.0 License. EEW, the spatial contiguity of the points with the variable values above their own threshold at the same time is a further request (as e.g. in Andreadis et al., 2005), which adds to the temporal contiguity typical of local extreme events definition (e.g. Hobday et al., 2016). 95

Characterisation
We introduce the characterisation of a EEWs based on two kinds of metrics: spatio-temporal indexes (sketched in Fig. 1) and strength indexes (sketched in Fig. 2).
The spatio-temporal indexes, useful to localise and describe an EEW in space and time (green shape in Fig.1), are: • the initiation, as the first day when at least one POT belonging to the EEW occurs; 100 • the duration T (yellow arrow in Fig. 1), as the time interval in which there are POTs included in the EEW. It is labeled by the maximum temporal difference between two POTs of the EEW, in day units; • the area A (grey area in Fig. 1), as the union of all the surface grid cells housing the POTs included in the EEW. It is labeled by the sum of these cells areas, measured in km2; • the width W, as the measure of the spatio-temporal region occupied by the EEW. It is computed as the sum over the 105 grid points covered by A of the spatio-temporal regions identified by the grid point area as base and the total time interval of POTs of the EEW referred to that grid point as height. It is measured in units of km2 × day; • the uniformity U, as the ratio between the width W and the spatio-temporal region defined by the prism with A as base and T as height:

=
(2.1) 110 It represents the percentage of the prism that is occupied by the EEW and quantifies how much (on average) the EEW is persistent on the single grid point belonging to A.
We excluded both EEWs with duration T<2 days (as e.g. Asch et al., 2019), to neglect possible transient spikes, and EEWs with area A < 4Δx × 4Δy (with Δx, Δy grid spacing in the zonal and meridional direction, respectively), since the estimated factor between effective resolution and grid spacing of the numerical models is 4 or more (Grasso, 2000). 115 The strenght indexes of a EEW can be defined starting from some quantities that are introduced locally (sketched in the top right box of Fig. 1). That is, considering the time series at each grid point (x,y), the j−th POT included in the EEW is characterised by the value Cj(x,y) of the ecosystem variable and by the intensity Ij(x,y) above the threshold p99(x,y) computed on the time series: ( , ) = ( , ) − 99( , ).
(2.2) 120 If we multiply Cj(x,y), Ij(x,y) and p99(x,y) by the local volume V(x,y), we obtain the local "masses" corresponding to the three quantities. Given the set J(x,y) of all the occurrence indexes j of the POTs that are referred to the specific grid point (x,y) and included in the EEW, we can define: https://doi.org/10.5194/bg-2020-34 Preprint. Discussion started: 27 February 2020 c Author(s) 2020. CC BY 4.0 License.
• the severity S, as the total mass of the variable supplied by the EEW, computed as the sum over the grid points covered by A of the local sum of the masses M(x,y) supplied by the POTs included in the EEW: 125 The severity is represented in a simplified way in Fig. 2a and is measured in kg.
• the excess E, as the total intensity above the "threshold" (i.e. the locus of points of the local thresholds, P99), associated to the EEW. Its formulation is analogous to Eq. (2.3), but referred to Ij (x,y) rather than Cj(x,y): Also the excess is represented in a simplified way in Fig. 2a and is measured in kg.
• the mean severity < S >, as the ratio between the severity and the width W of the EEW: Its unit is kg km-2 day-1; • the anomaly, as the ratio between the excess and the severity: 135 It represents the percentage of the excess in the severity of the EEW. This index, which is adimensional, is sketched in Fig. 2b for two different EEWs, which have same severity but different excess and, thus, anomaly. Since the locus of points of the thresholds of the second EEW is lower, this EEW has a higher value of the anomaly with respect to the first one and has a greater impact on the ecosystem. 140

Classification
The EEWs introduced in our formulation are identified starting from the local 99th percentile thresholds of the variable time series. The concept of "extreme" adopted in this work is related to the local characteristics of the marine ecosystem, which can be largely heterogeneous across the domain. Here we propose a classification of the EEWs, suitable to highlight the main features of the EEWs on the considered spatial domain. 145 For each index defined in Sect. 2.1.2, the values obtained for each EEW can be associated to all the points belonging to the A areas ( Fig. 1) and, then, a mean map on the spatial domain can be obtained by averaging point by point all the values of the related index. Finally, a fuzzy clustering analysis (Bezdek et al., 1984) can be conducted on the mean maps of all the indexes. In this way, different bio-regions of EEWs can be identified, depending on the relative weight of the indexes under consideration, and specific regimes of EEW can be therefore highlighted. 2015) following the online scheme described in Cossarini et al. (2017). The configuration in use has a horizontal resolution of 1/12°, with 75 vertical levels unevenly spaced. 155 In particular, we used the daily chlorophyll concentration computed at surface (i.e., averaged on the first 10 m). We restricted our investigation to the January-May period, since the occurrence of late winter-early spring blooms in the Mediterranean reproduces the heterogeneity of the marine ecosystem across the basin.

Chlorophyll EEWs: identification and characterisation
We applied the method illustrated in Sect. 2 to the surface chlorophyll concentration, so that C(x,y,t) ≡ chl(x,y,t), daily sampled, provided by the simulation of the Mediterranean Sea biogeochemistry in the 1994-2012 period (Sect. 2.2). 170 From the ecosystem point of view, the more suitable indexes to describe the phenomenology of the surface chlorophyll EEWs (i.e., exceptionally high and prolonged "blooms", as clarified in Introduction) are the mean severity <S>, the anomaly AN, the duration T and the uniformity U (Sect. 2.1.2). In fact, considering the chlorophyll as a proxy for the phytoplankton biomass (e.g. Boyce et al., 2010), the mean severity provides the mean amount of biomass supplied by the EEW to the surface layer in 1 day, over a unit area of 1 km2. The anomaly index instead represents the amount of chlorophyll 175 anomalously high with respect to the local history of the ecosystem. The duration T measures the overall ongoing impact on the marine ecosystem, where we intend "impacts on the ecosystem" as responses of ecosystem processes and status to an excess of phytoplankton biomass. On the other hand, the uniformity index quantifies the local persistence of the chlorophyll EEW on the points included in the area A. In fact, keeping constant the values of A and T, an EEW with higher U will affect the single unit of A for longer times, with higher potential ecological consequences on the ecosystem unit. 180 Considering the temporal extension of the simulation (approximately equal to 7000 days), the number of POTs in each grid point is by construction equal to 70. Mapping the 99th percentile threshold values computed at each grid point on the whole basin (Fig. 3), it can be noticed that grid points that are near in space exhibit small differences in their threshold values and also that different patterns are recognisable in the basin. Hereafter, we use the abbreviations indicated in Fig. 3  (1 km), it is able to capture a surface signature typical of deep convection dynamics (second column, compared with the first one, on 20th March). However, the comparison between model and satellite data points out how the cloud coverage affecting the remote sensing measurements is a limiting factor for the reconstruction of the spatio-temporal dynamics of the extremes of chlorophyll. Comparing the modelled chlorophyll maps (second column) with the patterns of the daily area A of the EEW (third column), we observe that the EEW patch actually includes points with noticeably high chlorophyll values in the 200 region. Nevertheless, A contains also points with chlorophyll values that are lower on the same absolute scale, yet higher than the local 99th percentile thresholds (as ensured by the procedure). Moreover, the EEW patches appear to be advected by the velocity field (third column) and to follow both the convection weakening (see plots in consecutive panels) and the patches of high nutrient concentrations in the previous days (by comparison with the right panel referred to the day before).
From Tab. 1, we quantify a mean severity equal to 1.389 kg km-2 day-1 and an anomaly index equal to 0.205% for this EEW. 205 In fact, it has been the most severe and the sixth most anomalous in the totality of the EEWs identified in the Mediterranean domain, as reproduced by our simulation. This states that the huge amount of chlorophyll that it supplied was also considerably higher than the local history of the impacted ecosystem. Moreover, even if the overall duration of the EEW was 17 days, each unit area was actually affected only for approximately 3 days (U T ≈ 3 days). This means that the EEW spread out in space and time with an articulated shape, as we can observe from Fig. 4. 210 In the Appendix, a further analysis conducted on three points that are located internally, externally and on the border of the EEW area showed that the EEW identification actually takes in account of all and only the relevant information associated to it and, thus, that the proposed method acts like a filter to properly circumscribe the extreme events in space and time. In fact, this specific EEW captured the dynamics of the exceptionally intense bloom observed in NWM in 2005 (Estrada et al., 2014;Mayot et al., 2016), triggered by a very strong vertical mixing (and deep convection, in the internal point), followed by the 215 restoring of stratification.
https://doi.org/10.5194/bg-2020-34 Preprint. Discussion started: 27 February 2020 c Author(s) 2020. CC BY 4.0 License.  Since most of the basin holds more than one EEW per year (Fig. 5a), the initiation time in each grid point and year was associated to the most severe EEW of that year. We found that the most severe chlorophyll EEW occurred mainly in the winter months in the central and southern open-sea part of the basin, later in the early spring period in NWM, central ALB, north ION, ADS, AEG and Rhodes Gyre (Fig. 5b). The duration of the chlorophyll EEWs covered up to 90 days in the southern part of the eastern basin, whereas it decreased to 30 days in NWM and Rhodes Gyre and to approximately 15 days 230 in ALB and ADS (Fig. 5c). High values of duration were typically associated to low uniformity (e.g., southern ION and LEV areas), while EEWs with high values of uniformity were found in ADS and ALB (Fig. 5d). The western Mediterranean displayed the EEWs with the highest mean severity, i.e. associated with the highest produced biomass, in ALB and NWM (Fig. 5e). Nevertheless, the regions with the highest anomaly are the eastern ADS and the northern Ionian Sea (Fig. 5f), despite their values of severity are around half of the ones of ALB or NWM. 235 Figure 6 displays the clusterisation provided by a fuzzy k-means analysis (Bezdek et al., 1984) conducted on the maps of the main indexes (i.e., duration, mean severity, uniformity and anomaly, Fig. 5), adopting a parameter of fuzziness equal to 2.
In Fig.6,  Along the simulated period, the western Mediterranean (except ALB), identified by clusters #3 and #7, did not display significant trends. On the other hand, in the Eastern sub-basin, (plus ALB), EEWs were longer and in ION and south eastern LEV (i.e. clusters #1, #2 and #4) also less uniform along the simulated years. A significant increase of anomaly (also higher than 0.2) was recorded in all the Eastern basin, except ADS and spotted areas in AEG and Rhodes Gyre (i.e. clusters #1, #2, #4, #5). 260

Sensitivity to the threshold
We conducted a sensitivity test of the method to different thresholds computed on the time series in each grid point. We repeated all the steps of the method (Sects. 2.1.1-2.1.3) for 98th and 99.5th percentile thresholds and we computed the mean value of the indexes in each cluster of Fig. 6. Finally, we computed the total means, averaging the mean values of the indexes on the seven clusters. Figure 7 shows the results, compared with 99th percentile (i.e., p99) reference threshold. 265 Duration and anomaly indexes show decreasing values for increasing thresholds. On the contrary, mean severity and uniformity display increasing values. The relative cluster ranks are generally preserved. Moreover, clusters #4, #6, #7 and #1 maintain the highest values of duration, uniformity, mean severity and anomaly, respectively, for all the selected thresholds, except in case of overall mean anomaly for 98th percentile, in which anomaly values of cluster #7 overcome the ones of cluster #1 (Fig. 9). 270

Discussion
In this work we propose a new method to tackle extreme events in the marine ecosystems on the basin scale. The method is then applied to the surface chlorophyll in Mediterranean open-sea areas to investigate maxima in the winter-spring blooms.
One of the key points of this method is the definition of "extreme event". In fact, the spatial extension of the extreme events is scarcely treated in literature, despite it can be an important ecosystem indicator, e.g. to predict a possible recovery of the 275 ecosystem (O'Neill, 1998;Thrush et al., 2005), and it is sometimes estimated a posteriori (Rabalais et al., 2002). On the contrary, the spatial contiguity of the local extremes has been here requested in addition to the temporal one, following Andreadis et al. (2005). In this way, the definition of the extreme from the time series in the grid point has been extended to https://doi.org/10.5194/bg-2020-34 Preprint. Discussion started: 27 February 2020 c Author(s) 2020. CC BY 4.0 License. define the Extreme Events Wave (EEW), which covers an extended area for a certain time duration. Consequently, the metrics necessary to characterise and classify the biogeochemical extremes, that can be introduced for time series of specific 280 sites (as e.g. in Hauri et al., 2013;Hobday et al., 2016;Asch et al., 2019;Salgado-Hernanz et al., 2019), have been then further developed to describe the shape and strength of the EEWs and to provide meaningful insights of the biogeochemical phenomenology.
In our specific application to the surface chlorophyll, the mean severity index associated to a chlorophyll EEW can be interpreted as the mean amount of biomass supplied daily to the sea surface over a unit area and could be used as an 285 indicator of eutrophication (Gohin et al., 2008;Ferreira et al., 2011) and of food availability for secondary production (Calbet and Agustí, 1999;Ware and Thomson, 2005). The map of the mean severity index obtained for the 1994-2012 period The anomaly feature, i.e. the case of supplied biomass much higher than usual for a certain area, can be instead ascribed to the inter-annual variability of the blooms (as in Mayot et al., 2016). As an example, the reconstruction of the most severe and On the other hand, the uniformity feature, i.e. the persistence of a chlorophyll EEW on a certain area, can be linked to 300 specific spatial constraints that circumscribe the EEW. In particular, the circulation structure can play an important role in providing the high values of uniformity of Fig. 5d. In fact, permanent cyclonic gyres in ADS (which impose also a topological constraint) and in northern LEV (i.e., Rhodes Gyre; Pinardi et al., 2015) potentially support a major vertical transport of nutrients and, consequently, higher values of biomass (Siokou-Frangou et al., 2010). Moreover, regular upwelling near the southern coast of Sicily can explain the high uniformity values in the Sicily Strait (e.g., Patti et al., 2010). 305 Finally, other spotted areas of high uniformity in ALB, SWW and TYR areas are characterised by semi-permanent mesoscale structures, associated to the inflow of Atlantic water (Navarro et al., 2011), to eddies originated from the Algerian Current (Morán et al., 2001) and to the dynamics of the northern TYR gyre (Artale et al., 1994;Marullo et al., 1994;Marchese et al., 2014), respectively.
The fuzzy k-means analysis used mean severity, anomaly, duration and uniformity to classify the EEWs in the Mediterranean Sea. The initiation index was excluded from the computation since in most part of the basin there is more than one EEW per grid point (Fig. 5a). However, as general characterisation of the EEWs occurrence, the initiation index shows a south-north gradient from winter to early-spring for the most severe EEWs in the Mediterranean open-sea (Fig. 5b) The clusterisation we obtained (Fig. 6) reveals a subdivision of the Mediterranean Sea that has several similarities to previous Mediterranean bio-regionalisations (D'Ortenzio and Ribera D'Alcalà, 2009;Lazzari et al., 2012;Ayata et al., 2018;Salon et al., 2019). This points out that the four indexes are meaningful in characterising the heterogeneity of the basin. 320 In particular, the cluster #7, corresponding roughly to the northwestern area (as "Bloom" region in D'Ortenzio and Ribera d' Alcalà, 2009;NWM in Lazzari et al., 2012), has been associated to the highest mean severity (Tab. 2). A decreasing gradient of the mean severity is observed toward the eastern Mediterranean areas (clusters #3 and #6 showed higher values of mean severity than clusters #1, #2, #4 and #5), in agreement with the west-to-east oligotrophication gradient (e.g., D'Ortenzio and Ribera d'Alcalà, 2009;Colella et al., 2016) and of the maxima of surface chlorophyll gradient (i.e. map of 325 amplitude index by Salgado-Hernanz et al., 2019). Moreover, this cluster is characterised also by a very high anomaly content, highlighting that blooms belonging to cluster #7 occasionally can supply a huge amount of chlorophyll, as in case of the already mentioned EEW in 2005, which occurred after a deep convection event (see Appendix). This interpretation is in agreement with the one referring to the "High Bloom" regime, which in some years takes the place of "Bloom" regime in NWM area (Mayot et al., 2016). 330 Cluster #1 identified the regime of the highest anomaly (i.e. of high inter-annual variability of the EEWs) and a decoupling between mean severity and anomaly. This regime in Mediterranean Sea is found in Northern ION and eastern ADS (Fig. 6).
Both cluster #5 and #6 are associated to high uniformity and low anomaly, i.e. to EEWs spatially well-localised and with regular occurrence over the years, respectively. Nevertheless, cluster #6 is characterised by EEWs with higher content of biomass (i.e. more severe) and running out more quickly (i.e. shorter) than #5. In this way, ALB, coastal SWW, ADS and 335 central part of Rhodes Gyre (i.e. cluster #6, Fig. 6) are differentiated by south western ION, AEG, outer part of Rhodes Gyre (i.e. cluster #5).
Cluster #4 displays the longest and less uniform chlorophyll EEWs, which are few (not shown), consistently with the relative high anomaly, and have low severity. This typology of EEWs identifies a regime of spatially diffuse blooms whose chlorophyll values do not differ markedly from the chlorophyll means in the concerned area. This cluster covers a large part 340 of south eastern Mediterranean Sea (Fig. 6), crossed by the Atlantic-Ionian Stream and the Cretan Passage Southern Current . We ascribe the very low uniformity and the high overall duration of these EEWs to the transport and spreading of chlorophyll along the meanders of these currents (not shown).
Finally, clusters #2 and #3 display in-between conditions with respect to others, since the values of all the indexes are intermediate, except the mean severity, which is very low in cluster #2 and relatively high in cluster #3. In the Mediterranean 345 https://doi.org/10.5194/bg-2020-34 Preprint. Discussion started: 27 February 2020 c Author(s) 2020. CC BY 4.0 License. basin, this corresponds to the decreasing gradient in the severity between the central part of the western (i.e. cluster #3) and of the eastern (cluster #2) Mediterranean Sea (Fig. 6).
The obtained clusterisation (Fig. 6) has been used also to evaluate the long-term evolution of the ecosystem phenomenology.
Our results did not show any increase in the intensity (i.e. in the severity index) of surface chlorophyll EEWs in any clusters over the period 1994-2012 (Tab. 2). This result is in agreement with estimations of trends in the amplitude index by Salgado-350 Hernanz et al. (2019), except in NWM. Moreover, no significant trend (here defined as an annual variation of an index which was higher than 1%) has been estimated for none of the four indexes in the central and north western sub-basin. On the contrary, the eastern Mediterranean (plus ALB) showed trends in duration, uniformity and anomaly of chlorophyll EEWs. In particular, positive trends of duration found in areas with very uniform EEWs suggest a persistence of the blooms which has been prolonged in time. On the other hand, positive trends of duration, along with low values of uniformity (with also 355 negative trends) and higher anomaly, denote an increase of EEWs with articulated shapes in lower productive areas. The trends recognised in the eastern Mediterranean Sea suggest a possible higher tendency of this sub-basin to changes in the identified regimes, despite the productivity is lower than in the western Mediterranean. This is one of the features that could emerge only accounting for the local thresholds in the identification of the EEWs (Sect. 2.1.1).
We have applied the method to the surface chlorophyll, as one of the more representative and investigated variables of the 360 marine ecosystem, identifying maxima in the blooms as events which potentially influence the ecosystem function (e.g., food web and carbon fluxes). Other variables whose impacts on the ecosystem can be relevant, such as HAB-like phytoplankton groups (Vila and Masó, 2005), can be investigated with our method, provided the availability of a continuous dataset in time and space. Moreover, the formulation illustrated in Sect. 2.1 could be extended to the full 4D case (i.e. to variables C(x,y,z,t)), adding the vertical dimension in the definitions of the local extreme (i.e., the POTs could be defined in each point 365 in a 3D space) and of the EEW (i.e., as 3D spatial volumes connected in time). The spatio-temporal indexes would refer to the spatial volume, instead of the area A, and to 4D width and prism associated to the definition of uniformity. The 4D formulation could be applied to investigate the marine hypoxia, identifying volumes in time with extremely low values of oxygen. In this case, a proper threshold for the local extremes would be fixed as the same value at each point, in connection with the impacts on benthic fauna and fish species, which have physiological limits (e.g., Rabalais et al., 2002, Vaquer-370 Sunyer andDuarte, 2008). Strength indexes would be modified and defined as vertical profiles, instead of scalar metrics, to show the intensity and the depth of the bottom ecosystem stress. Therefore, the novelty of our method (i.e. the temporal and spatial connection of extremes) allows to compute the extension of the (connected) spatial 3D volumes under hypoxia to estimate the fish survival probability, enhanced by swimming (avoidance) behaviour (Rose et al., 2017).
A critical parameter of our method is the choice of the local percentile threshold (Sect. 2.1.1). In our case, it was computed 375 as the 99th percentile on the surface chlorophyll time series in each grid point. A priori, the choice of a higher (lower) threshold corresponds to a definition of "extreme" value which is narrower (broader) than the reference one. As shown in Fig. 7, the choice of higher thresholds increases the mean severity index, since it is computed on local values which were higher. On the contrary, both the anomaly and the duration indexes decrease at higher thresholds, because of the occurrence https://doi.org/10.5194/bg-2020-34 Preprint. Discussion started: 27 February 2020 c Author(s) 2020. CC BY 4.0 License. of local POTs for a smaller number of days (i.e., lower duration) and of a less effective detectability of the inter-annual 380 variability (i.e., lower anomaly). The increase of the uniformity index is due to the promotion of grid points with more similar values of higher local threshold (i.e. closer in space, see Fig. 3). However, uniformity shows lower sensitivity to the threshold than the other indexes because of occurrence of POTs in few grid points with thin spatial connectivities extending for great distances (as discussed in Sheffield et al., 2009). In this case, a further sensitivity analysis on the area covered by the EEWs could be envisaged to identify a minimum area threshold (stricter than the 4Δx × 4Δy constraint introduced in 385 Sect. 2.1.2) to better characterise the uniformity property.
Overall, Fig. 7 shows that the identification of the clusters with the highest content of all the indexes has been generally maintained both in case of higher and lower thresholds, proving that the main regimes of chlorophyll EEW (corresponding to clusters #1, #4, #6, #7 of Fig. 6) have been individuated in a robust way. Only the intermediate regimes (clusters #2, #3, #5) could display other relative weights of the indexes for threshold higher/lower than the reference. 390 Since different variables of interest could highlight a different sensitivity of the indexes, we think that carrying on analyses with different local thresholds could help to identify the specificities of the phenomenology underlying the extremes.

Conclusions
The present study provides a methodology to describe the statistical extremes in the marine basin-scale ecosystem, supported by an ecological interpretation. 395 A key issue of the method is the request of contiguity both in time and in space of the peaks over the local threshold of the ecosystem variable. This constraint allowed to define individual events as Extreme Events Waves (EEWs) occurring in localised spatio-temporal regions. In particular, we accounted for the contiguity of the local extremes, which is an aspect that has been rarely considered in literature. At the same time, our choice to start from local thresholds, computed as a percentile of the time series in the grid point, allowed to maintain a definition of "extreme" relative to the local ecosystem properties. 400 For a biogeochemical variable evolving in a two-dimensional space, we proposed a set of indexes for EEWs to describe their initiation, duration, total covered area and (spatio-temporal) uniformity, as well as their (mean) severity and anomaly, as measures of overall intensity and inter-annual variability, respectively.
In the specific application to the open-sea Mediterranean chlorophyll, we characterised the maxima of "blooms" (in a local and statistical sense), as EEWs which potentially influence the ecosystem function. A cluster analysis conducted on the 405 indexes associated to the covered areas allowed to identify four main regimes. We associated the occurrence of blooms with high mean severity and high inter-annual variability to the north western Mediterranean Sea; blooms with high inter-annual variability (associated to intermediate intensity) to the northern Ionian Sea; regular and spatially well-localised blooms in Alboran and south western Mediterranean Sea, south Adriatic Sea and Rhodes Gyre; weak and diffuse blooms in the south eastern Mediterranean Sea. 410 https://doi.org/10.5194/bg-2020-34 Preprint. Discussion started: 27 February 2020 c Author(s) 2020. CC BY 4.0 License.
We did not observe significant trends (i.e. annual variations higher than 1%) of the mean severity of chlorophyll EEWs across the Mediterranean basin, despite some trends are found for other indexes.
Comparison of the results with available data and previous studies supports the reliability of the method, which could be promisingly applied to other ecosystem variables. However, sensitivity analyses are recommended to select suitable thresholds to highlight the typology of the extremes under consideration. In the internal point A, 30-40 days preceding the EEW onset were characterised by strong heat losses up to 1000 W/m2 (top panel of Fig. A.2) due to the wind field (not shown). This led to a strong deep convection that mixed all the water column, down to the sea bottom (second panel of the same figure). Surface and subsurface nitrate, whose concentration at the beginning of the year was already above the climatological values, was further enhanced during the mixing (fourth panel). 425 As soon as the stratification was quickly established (second panel) and the surface temperature rose (third panel), an abrupt rise of surface chlorophyll occurred (bottom panel). The surface chlorophyll, that in January, February and the first half of March had exhibited values much lower than the climatological ones due to the strong convective phase, in the third week of March increased by a factor of almost 800% in 4 days. A full consumption of surface nitrate can be observed in the same days (fourth panel). A following weaker mixing phase (in the half of April) replenished the surface layers with a relatively 430 low amount of nitrate (yet above their climatological values), triggering two weak episodes of increasing chlorophyll. On the whole, the features here described are in agreement with the characterisation of the chlorophyll blooms in the NWM area and, in particular, of the 2005 event (Barale et al., 2008;Estrada et al., 2014;Mayot et al., 2016).
The interpretation of the results referred to the peripheral point B belonging to the EEW area (Fig. A.3) is similar to the previous considerations about the internal point A, but in presence of a less intense vertical mixing. 435 On the other hand, in the point C (Fig. A.4), external to the EEW area, an evident stratification of the water column below the 30 m of depth has been maintained for all the winter months (January-February-March), despite the cooling of the surface layers. The nitrate content in the surface and subsurface layers was much lower with respect to the deeper ones and only a small increase of the surface chlorophyll developed during the duration of the EEW.
Therefore, the strong deep convection (related to the inter-annual variability of the local vertical mixing) appears to be the 440 key factor for the exceptionality of this EEW. It is worth to remind that our method identified the spatio-temporal region https://doi.org/10.5194/bg-2020-34 Preprint. Discussion started: 27 February 2020 c Author(s) 2020. CC BY 4.0 License. The excess (dark green portion, on the right) is the part of this volume that is above the locus of points of the 99th percentile threshold, i.e. P99, delimited by the orange contour. Fig. 2b): The anomaly, as the percentage of the excess with respect to the severity, is compared between two cases, which are referred to two EEWs with the same severity, but with different P99 locuses of points.       Table 2: Mean and standard deviation of duration, uniformity, mean severity and anomaly indexes within the seven clusters in Fig. 6. Pale red (blue) colour refers to an annual increase (decrease) higher than 1% in the 1994-2012 period.