Spatial and seasonal variability of horizontal temperature fronts in the Mozambique Channel for both epipelagic and mesopelagic realms

Introduction Ocean fronts are moving ephemeral biological hotspots forming at the interface of cooler and warmer waters. In the open ocean, this is where marine organisms, ranging from plankton to mesopelagic fish up to megafauna, gather and where most fishing activities concentrate. Fronts are critical ecosystems so that understanding their spatio-temporal variability is essential not only for conservation goals but also to ensure sustainable fisheries. The Mozambique Channel (MC) is an ideal laboratory to study ocean front variability due to its energetic flow at sub-to-mesoscales, its high biodiversity and the currently debated conservation initiatives. Meanwhile, fronts detection relying solely on remotely-sensed Sea Surface Temperature (SST) cannot access aspects of the subsurface frontal activity that may be relevant for understanding ecosystem dynamics. Method In this study, we used the Belkin and O’Reilly Algorithm on remotely-sensed SST and hindcasts of a high-resolution nested ocean model to investigate the spatial and seasonal variability of temperature fronts at different depths in the MC. Results We find that the seasonally varying spatial patterns of frontal activity can be interpreted as resulting from main features of the mean circulation in the MC region. In particular, horizontally, temperature fronts are intense and frequent along continental shelves, in islands’ wakes, at the edge of eddies, and in the pathways of both North-East Madagascar Current (NEMC) and South-East Madagascar Current (SEMC). In austral summer, thermal fronts in the MC are mainly associated with the Angoche upwelling and seasonal variability of the Mozambique current. In austral winter, thermal fronts in the MC are more intense when the NEMC and the Seychelles-Chagos and South Madagascar upwelling cells intensify. Vertically, the intensity of temperature fronts peaks in the vicinity of the mean thermocline. Discussion Considering the marked seasonality of frontal activity evidenced here and the dynamical connections of the MC circulation with equatorial variability, our study calls for addressing longer timescales of variability to investigate how ocean ecosystem/front interactions will evolve with climate change.


Introduction
Ocean fronts are ubiquitous discrete features in the ocean which are now more readily available thanks to the everincreasing spatial resolution of satellite data Kartushinsky and Sidorenko, 2013;Belkin, 2021). Fronts are defined as a sharp transition between waters of different properties (temperature, salinity, chlorophyll-a, nutrients, etc.) and are usually detected by high local horizontal gradients in those same properties (Largier, 1993;Cayula and Cornillon, 1995;Shimada et al., 2005). Frontal phenomena encompass a wide range of space and time scales, from submesoscale structures to thousands of kilometers and from hours to hundreds of years (McWilliams, 2021;Belkin, 2021). More specifically, some fronts are seasonally persistent (Acha et al., 2004;Vinayachandran et al., 2021;McWilliams, 2021), while others display clear seasonal patterns (Park and Chu, 2006;Castelao et al., 2006;Bosse and Fer, 2019).
Indeed, submesoscale fronts were previously thought to be confined in the mixed layer (Fox-Kemper et al., 2008) but recent studies show that fronts could penetrate much deeper than the mixed layer, and reach depths up to 1000m (Bettencourt et al., 2012;Yu et al., 2019;Siegelman et al., 2020b). There is thus an opportunity to infer deep circulation as well as verticallyintegrated biological hotspots from surface mesoscale activity. This is essential for understanding the spatio-temporal variability of these critical ecosystems, at both ocean surface and interior, in order to design conservation management strategies and ensure sustainable resources (Alemany et al., 2009;Nieblas et al., 2014;Watson et al., 2018;Scales et al., 2018).
In this study, we investigate the temperature fronts in the Mozambique Channel (MC). The MC is part of the greater Agulhas system (Lutjeharms, 2006) starting North of Madagascar and reaching south all the way to South Africa ( Figure 1). The circulation in the Mozambique Channel is highly energetic and unstable. The South Equatorial Current (SEC) parts in two at around 17°S east of Madagascar forming the Northeast and Southeast Madagascar currents (SEMC and NEMC, respectively). The MC is characterized by a highly variable Mozambique Current flowing southward along the Mozambican shelf (DiMarco et al., 2002;Lutjeharms et al., 2012) and a succession of mesoscale eddies propagating southward at a speed of about 3 to 6 km/day (de Ruijter et al., 2002;Halo et al., 2014). Eddies in the MC have a typical scale of several hundreds of kilometers (de Ruijter et al., 2002;Halo et al., 2014;Ternon et al., 2014), and a frequency of occurrence of 4 to 7 per year (Schouten et al., 2003). Those eddies display a strong barotropic component and some are felt all the way to the bottom of the channel (i.e., over -2000m) (de Ruijter et al., 2002;Ridderinkhof and de Ruijter, 2003;Schouten et al., 2003;Halo et al., 2014). The most intense eddies are anticyclonic eddies (AE), also referred to as Mozambique Rings. They form at the narrowest part of the MC around 12-16°S and travel southward along the Mozambican shelf (Halo et al., 2014), where the flow is most intense (de Ruijter et al., 2002) and where eddy kinetic energy (EKE) is the highest (Biastoch and Krauss, 1999). Cyclonic eddies (CE) are more numerous, yet smaller than AE in the channel (Halo et al., 2014). They usually form at the edge of bigger AE on the eastern side of the channel and south of 20°S, to then travel to the southwest (Halo et al., 2014;Roberts et al., 2014). The formation of submesoscale filaments and fronts at the edge of mesoscale eddies is a common feature in turbulent regions and has been extensively documented from satellite imagery (Munk et al., 2000;Zhang and Qiu, 2020) and from observations (Read et al., 2007;Roberts et al., 2014;Adams et al., 2017;Siegelman et al., 2020a;Tarry et al., 2021). The highly energetic flow at sub-to-mesoscale of the MC makes it thus ideal to study fronts.
The MC also exhibits a strong seasonal cycle and experiences biannual monsoons north of 20°S (Schott et al., 2009). During the northeast monsoon, from December to March, winds blow from the northwest, while winds blow from the southeast during the southwest monsoon, from June to November. South of 20°S, winds blow from the southeast all year round (Obura et al., 2019). The monsoon reversals of the north MC also modulate currents intensity (Ridderinkhof et al., 2010;Obura et al., 2019) and create a distinct seasonal pattern, which affects the large anticyclonic eddies forming at 12-16°S (Donguy and Piton, 1991). Finally, altimetry-derived EKE climatology for different parts of the MC ( Figure 1) shows a strong seasonal variability, especially north of 25°S. The MC energetic turbulent flow, seasonal cycle and exceptional biodiversity (Obura et al., 2019) thus motivate the investigation of temperature front seasonal variability.
Historically, fronts have mostly been detected from satellite imagery. Originally, fronts were simply identified manually from SST images (Legeckis, 1978), but more sophisticated automatic front detection algorithms have been developed since then. Front detection algorithms are usually based on SST gradients, even though some also use entropy (Shimada et al., 2005), neural networks (Torres Arriaza et al., 2003) or singularity exponents (Turiel and del Pozo, 2002;Turiel et al., 2008;Tamim et al., 2015). The Canny algorithm (Canny, 1986) is often referred to as the original "gradient method" and is based on the high gradient property of fronts. However, this method often results in enhanced noise from the differentiation (Belkin and O'Reilly, 2009). The Cayulla-Cornillon Algorithm (CCA) addresses this shortcoming by operating with histograms to detect the edge between two water masses (Cayula and Cornillon, 1995). This method requires fine-tuning to find the appropriate threshold distinguishing fronts from background noise, but has been successfully used in many different systems: in the Baltic Sea (Kahru et al., 1995), in the Atlantic Ocean (Mantas et al., 2019) or in the South China Sea (Wang et al., 2001). The performance of the CCA was further improved by adding multiple sliding windows (Nieto et al., 2012) and then applied in the Indian (Nieblas et al., 2014) and Pacific (Nieto et al., 2017) Oceans. As opposed to histogram-based methods, the algorithm proposed by Belkin and O'Reilly is based on SST gradients and effectively remove the noise generally produced by other gradient methods, while preserving the shape of the detected front (Belkin and O'Reilly, 2009). The Belkin and O'Reilly algorithm (BOA) produced fronts dataset in the South China Sea (Wang et al., 2001), East China Sea (Lin et al., 2019) and the Kuroshio current (Liu et al., 2018). Fronts can originate from the interaction between strong currents and topography anomalies, such as shelf-breaks, islands and seamounts (Heywood et al., 1990;Condie, 1993;McWilliams, 2021); at the edge of eddies in turbulent areas (Munk et al., 2000;Read et al., 2007;Siegelman et al., 2020a) and in upwelling zones (Bettencourt et al., 2012;Rossi et al., 2013;Sabarros et al., 2014;Calil et al., 2021). The Island Mass Effect (Doty and Oguri, 1956) in particular plays an important role in frontogenesis, as it has been shown to be associated with the generation of vortices in the lee of islands (De Falco et al., 2022). Although it is very rare to attribute clearly a single driver to frontogenesis, it is instead often a complex combination of favorable conditions (McWilliams, 2021).
The complex dynamic of processes behind frontogenesis (McWilliams, 2021) suggests that fronts detection relying solely on remotely sensed SST might miss important front features and variability at depth, and estimates of the magnitude of SST gradients in satellite products is highly sensitive to resolution and interpolation technique (Vazquez-Cuervo et al., 2013). In this study, we use remotely-sensed SST and outputs of a highresolution realistic ocean model to investigate the spatial and seasonal variability of thermal fronts at different depths in the Mozambique Channel. The paper is organized as follows: we first introduce our material and methods for front detection and we follow with a description of front variability in the Mozambique Channel. The article ends by a discussion and a few concluding remark 2 Material and methods

Domain subregions
Before studying thermal front variability, we divided our domain into 4 subregions, as seen in Figure 1. The rationale for this meridional subdivision is based on a literature review about drivers of oceanic variability (i.e. bathymetry effects, monsoon forcing) in the channel and relies on the dependence of eddy characteristics and therefore fronts to the Coriolis force. The "Islands" subregion (10-13.5°S, 38-49.5°E) starts just above the Comoros archipelago and extends to 13.5°S at the entrance of the Angoche Basin. The "Northern Mozambique Channel" (NMC) subregion (13.5-20°S, 34-48°E) ends at 20°S where the influence of monsoons weakens substantially (Schott et al., 2009). The "South Mozambique Channel" (SMC) subregion (20-25°S, 33-45°E) goes as far south as 25°S, which corresponds to the southern tip of Madagascar island. Finally, the "Agulhas" subregion (25-30°S, 31-50°E) covers the southernmost part of our domain from 25°S to 30°S.

Model configuration
In this study, temperature fronts are derived from a ROMS-CROCO (Regional Ocean Modeling System and Coastal and Regional Ocean COmmunity model) configuration of the MC (Shchepetkin and McWilliams, 2005). ROMS-CROCO solves the primitive hydrostatic equations for momentum and tracer variables and successfully resolves regional and coastal oceanic processes at high resolution, as well as mesoscale eddies (Halo et al., 2014). The model configuration consists in 3 nested grids (two-way nesting) with the child domain at the highest resolution of 1/36°used here. 60 vertical sigma (topography following) layers are used. The bathymetry used in the model is from GEBCO 2014 (Weatherall et al., 2015). A bulk formulation (Fairall et al., 1996) is used to compute surface turbulent heat to force the model, based on ERA-Interim Reanalysis daily mean atmospheric fields (2 m air temperature, relative humidity, surface wind speed, net shortwave and downwelling longwave fluxes, and precipitations). Wind stress forcing is also from the ERA-Interim Reanalysis. At the lateral open boundaries, mixed active-passive conditions are used (Marchesiello et al., 2001) with forcing data from GLORYS reanalysis at 1/4°resolution. The generated daily output covers 22 years, from 1993 to 2014. This 1/36°ROMS-CROCO configuration of the Mozambique Channel, referred to as MOZ36, was first used by Miramontes et al. (Miramontes et al., 2019) to document the deep oceanic circulation in the Zambezi turbidite system (SMC subregion in Figure 1)

Climatological and satellite datasets
We used climatological seawater properties from CARS2009 (Ridgway et al., 2002) generated from compiled in-situ observations spanning from January 1941 to May 2009, on a 1/ 2°grid. We also retrieved World Ocean Atlas 2018 annual mean temperature  and salinity . Remotely-sensed Sea Surface Height (SSH) was taken from Global Total Surface and 15 m Current (Rio et al., 2014, COPERNICUS-GLOBCURRENT) (Chin et al., 2017). In this work, we use the full MUR dataset from January 2003 to December 2020, for the MC. SST MUR data comes with 2 data quality flags. Flag 1 corresponds to the estimated error standard deviation of the SST data and flag 2 is the time lag to most recent 1 km data. Flag 1 is available for the entire dataset while flag 2 is only included after July 2016.

Numerical methods
We use the Brunt-Väisälä frequency N as a proxy for stratification. N is given by: Where g is the gravitational constant, r 0 is a reference density set at 1025 kg.m -3 r is the seawater density and z its depth.
Intra-seasonal and inter-annual eddy kinetic energy (EKE) were computed following a similar method to Dewitte et al. (Dewitte et al., 2021). First, intra-seasonal current anomalies (u IS and v IS ) were estimated by subtracting monthly mean currents to the daily velocity dataset. Then, the mean intra-seasonal EKE was computed as: Second, inter-annual current anomalies (u IA and v IA ) were computed by subtracting the climatology to the daily velocity dataset. Then, the mean inter-annual EKE was computed as: Following the same methodology, MOZ36 intra-seasonal EKE was also computed at -100 m, -500 m and -1000 m.
We chose the Belkin and O'Reilly algorithm (Belkin and O'Reilly, 2009) to detect the location and intensity of temperature fronts in the MC. More specifically, we use Ben Galuardi's BOA pseudo-code written in R and publicly available on https://github.com/galuardi/boaR (accessed on 19 May 2022). BOA first applies a contextual median filter until convergence to reduce noise originating from natural variability and artifacts. The algorithm then uses 3x3 kernel Sobel operator to compute the gradient magnitude (Belkin and O'Reilly, 2009). BOA temperature gradients were computed on MUR SST but also on MOZ36 temperatures between 0 and -1000 m, leveraging the fact that the model is not limited to the surface, as would be satellite datasets. From the surface to -300 m, the model outputs were interpolated on a 10 m vertical resolution grid, while below -300 m we use a resolution of 50 m.
We derive several statistics to characterize temperature fronts and to evaluate their variability. BOA being a gradient-based method of front detection, the algorithm produces a continuous field of temperature gradient magnitude, or front intensity in°C/ km. As such, a threshold must be determined to define what is considered a front. Choosing a high threshold might prevent weak fronts from being recorded, while a weak threshold might blur any meaningful signal into noise. We thus define front pixels as pixels whose temperature gradient magnitude is equal or above the chosen threshold. We finally define the probability of front occurrence per month for each ocean pixel as the percentage of days when the gradient value qualifies as "front", i.e. when the pixel gradient magnitude is over the threshold, over the total number of days for each month of simulation. Finally, we computed sea level anomalies (SLA) from the mean MOZ36 sea-surface heights in order to link frontal activity with (sub)-mesoscale features in the MC. And dominant timefrequencies were studied based on 1-D discrete fast Fourier transforms (FFT) computed with the scipy.fft package.

Model validation 2.5.1 EKE
To validate MOZ36 (sub-)mesoscale structures on intraseasonal and inter-annual time scales, we evaluated the mean eddy kinetic energy (EKE) derived from the model and remotely-sensed currents (see supplementary material Figure  S1). The spatial EKE patterns are similar in both model and satellite data: the EKE is most intense on the western part of the MC and in its southern part, in the wake of the SEMC. Note that, the slightly higher EKE from the model could be well-explained by its higher resolution.

Hydrography
When comparing the SST variability in the MC from MOZ36 and from MUR SST ( Figure S2), both datasets exhibit comparable seasonal and spatial SST patterns. Additionally, and to further validate seawater properties of MOZ36, we compared World Ocean Atlas 2018 temperature  and salinity  annually-averaged vertical profiles on a meridional section at 42.125°E, in the center of the channel ( Figure S3). Observed and modelled vertical profiles show good consistency, confirming that MOZ36 is suitable to study temperature fronts in the MC.

Stratification
We used CARS2009 (Ridgway et al., 2002) temperature and salinity profiles to evaluate the realism of the modelled vertical stratification thanks to the Brunt-Väisälä frequency from the surface to -1000 m. The comparison between the mean model and CARS2009 squared Brunt-Väisälä frequencies for the 4 MC subregions ( Figure S4) reveals that MOZ36 is able to reproduce a range of values similar to observations, with a maximum N 2 of about 0.2*10 -3 s -2 . For both datasets, N 2 in the Agulhas region is markedly smaller than the other 3 subregions and peaks at -60 m. The Islands, NMC and SMC subregions peaks between -60 and -100 m in the model and between -80 and -125 m in CARS2009. Overall the model reproduces a rather realistic mean thermocline structure, providing high confidence in its ability to simulate well the main aspects of ocean dynamics.

Front threshold definition
We define a single temperature front threshold for all depths using statistical distribution and quantiles of MOZ36 temperature front intensity for all ocean pixels across the 22 years of simulation and 4 depth levels: surface, -100 m, -500 m, and -1000 m ( Figure  S5). Surface temperature gradients go up to 2.3°C/km, and 70% of those are below 0.05°C/km. We subsequently define front pixels as pixels whose temperature gradient magnitude is equal or above 0.05°C/km, and the probability of front occurrence per month for each ocean pixel as the percentage of days when the pixel gradient magnitude is over 0.05°C/km, over the total number of days for each month of simulation.

Spatial variability 3.2.1 Surface
The mean surface temperature front intensity and probability maps (Figures 2A, B) reveals that SST fronts derived from the model are more intense and more likely to appear along the continental shelves, in the wake of both Northeast and Southeast Madagascar currents (NEMC and SEMC) and more generally southward of 25°S.
The equivalent mean maps of temperature front intensity and probability derived from the MUR dataset exhibit a spatial distribution similar to that of the modelled fronts ( Figures 2C,  D), despite a slightly lower range of values. Figure 3 represents the mean vertical profiles of front intensity in°C/km from MOZ36, the shading around the mean corresponds to the seasonal standard variation. The vertical profiles reveal three distinct regimes, similar in all 4 subregions of the MC:

Ocean interior
1. From the surface to about -100 m, the frontal activity increases abruptly to reach its maximum to 0.75-0.9°C/ km at -100 to -200 m, in the vicinity of the mean thermocline. 2. From -200 m to -500 m, the front intensity decreases sharply to reach 0.35°C/km in the NMC region and around 0.05°C/km in the Islands, NMC and Agulhas regions. 3. From -500 m and deeper, frontal intensities decrease at a slower rate with depth. At -1000 m the NMC region exhibits the lowest frontal activity, followed by the Islands regions. The SMC and the Agulhas region show less of a vertical difference between surface and depth, than the NMC and the Islands regions.
The vertical distribution of temperature front intensity is very similar to MOZ36 and CARS2009 stratification profiles ( Figure S4). Indeed, the Islands (Agulhas) subregion has the A B D C  Seasonal variability ceases to be significant, i.e. it amounts to less than 5% of the front intensity mean, below -150 m for the Islands region, below -350 m and -450 m for the NMC and SMC regions (respectively) and below -70 m in the Agulhas region.
Frontal activity peaks at around -100 m, in the vicinity of the mean thermocline ( Figure 3). The resulting 22-year map of mean front intensity at -100 m ( Figure 4C) indicates stronger frontal activity than at the surface, as expected from the front intensity distribution ( Figure S5). The frontal activity is specifically stronger along the continental shelves, in the wake of the NEMC and SEMC and in the central part of the MC, in between Bassas da India and Juan de Nova islands across the Angoche Basin.
-500 m and -1000 m temperature fronts are most intense along the Mozambican shelf, near Maputo and in the path of the NEMC and SEMC (Figure 4).

Front typology
Selected snapshots of temperature front intensity at different depth and time of the year allowed us to distinguish three types of temperature fronts based on their signature behaviors ( Figure 5).

• Topographic fronts, as defined by McWilliams
McWilliams, (2021) are formed where currents encounter land masses or other sharp topographic changes such as seamounts. In this study, topographic fronts are found along Mozambique and Madagascar shelves and around islands (French Scattered Islands, Comoros Archipelago), they reach maximum intensity at -100 m. NEMC and SEMC paths can clearly be identified from the front intensity patterns at the north and south tips of Madagascar Island (Figure 4). • The next two types of fronts correspond to the "deformation fronts" in McWilliams's typology. They form at the edge of eddies.
• Anticyclonic deformation fronts (ADF) are located at the edge of anticyclonic eddies (AE) in the MC. AE usually form at the entrance of the Angoche Valley at 15°S and travel southward close to the Mozambican shelf to then merge or dissipate around 22°S and feed the Agulhas Current. ADF reach their maximum intensity at around -100 m and then weaken substantially at -500 m, to disappear completely by -1000 m.
• Cyclonic deformation fronts (CDF) form at the edge of cyclonic eddies (CE), generally south of 23°S. Most of those CE are driven by the SEMC: some go slightly northward up the southern tip of the Madagascar shelf, but most go southwest following the SEMC track towards the Agulhas current. CDF reach their maximum intensity at -500 m and can still be detected at -1000 m.
A summary of the front typology and locations can be found in Table 1.

Seasonal variability 3.3.1 Islands and NMC subregions
In the Islands (10°S-13°S) and NMC (13°S-20°S) regions, SST fronts are strongest in the wake of the NEMC between the Mean temperature front intensity for each MOZ36 subregion in°C/km. The shading shows the seasonal standard deviation from the mean. The dashed horizontal lines mark the depth below which the seasonal standard deviation of temperature front intensity is lower than 5% of its mean.  Note that the surface climatology of front intensity from MUR SST (Figure 7) shows little seasonal variability compared to the model. This aspect is explored later in the discussion part, where quality flags explain well why the seasonal cycle is significantly dampened compared to the model.
In the Islands subregion, the -100 m front intensity is greater than at the surface and peaks in late austral autumn (June). The signature 2-peak of the surface signal seems to have disappeared entirely, leaving only one clear maximum in winter. In the NMC area, the -100 m front intensity is higher than at the surface and peaks in austral summer, with a second local maximum in austral winter. At -500 m and -1000 m, there is not much seasonal variation anymore (Figure 7).

SMC and Agulhas subregions
In the SMC (20°S-25°S) and in the Agulhas (25°S-30°S) subregions, surface temperature fronts are more intense along the continental shelves and in the wake of the SEMC. The FFT analysis of the 22-year daily time series of SST front intensity from MOZ36 in the two southernmost regions reveals both 6month and 1-year periods, but this time the yearly signal exhibits stronger magnitude than the biannual one ( Figure 6). The climatology of surface temperature front intensity (Figure 7) shows stronger frontal activity during austral winter and spring, with a slight increase in December in the SMC.
Similarly to the NMC, the -100 m temperature front climatology in the SMC is stronger than at the surface and reveals two local maxima with comparable range, in austral summer and winter (Figure 7). In the Agulhas region, the -100 m front climatology (Figure 7) shows a much weaker seasonal variation, also seen on the FFT analysis ( Figure 6), with a slight maximum in austral winter.
At -500 m and -1000 m, little seasonal variation is seen in frontal activity (Figure 7). However, even though the fronts become less intense with depth, -500 m and -1000 m fronts in the SMC and Agulhas subregions are still more intense than in the Islands and NMC regions. Indeed, the interaction between the SEMC and local topography fosters the formation of cyclonic eddies and thus of CDF, which have shown to have a strong deep signature in the model ( Figure 5).

Data and method
The extended MOZ36 validation process conducted in this work and in Miramontes et al. (2019) has confirmed its suitability in depicting realistic basin-scale and regional dynamics of the MC. Yet, the temperature front seasonal variability exhibited by MUR SST data does not match perfectly MOZ36 results, despite consistent spatial patterns. This disparity can be explained by MUR SST cloud cover adjustments. A preliminary comparative study between MUR SST flags values and corresponding BOA temperature front intensity revealed that periods with higher cloud cover (mostly during austral summer) resulted in dampened frontal activity: the frontal structures are less sharp or non-existent, and the gradient magnitudes much less intense. These higher cloud cover areas correspond to flag 1 values being above 0.39°C and flag 2 values being above 20 hours. For example, Figure 8 shows snapshots of front intensity (C, D) and corresponding flag 1 (E, F) and flag 2 (G, H) for the first day of February and July 2020. One can distinctively see a difference in Figure 8C between the Mozambique shelf north of 17°S where fronts are dimmed, and the center of the MC at around 20°S, where the structures are much sharper. The top panels of Figure 8 show the climatology of monthly percentage of ocean pixels having acceptable flag values to study frontal activity, based on flag 1 for 2003-2020 ( Figure 8A) and on both flags for 2017-2020 ( Figure 8B). If we only take flag 1 into account, 50-90% of the dataset is acceptable, while if we consider both flags then only 10-50% of the dataset is suitable to study seasonal variability of submesoscale fronts in the MC. Regardless of how many flags are considered, a bigger proportion of pixels are not suitable in austral summer, compared to austral winter, providing a good explanation as to why so little seasonality is seen from MUR SST front intensity in Figure 7. Indeed, the difference between seasons in frontal activity could be significantly dampened by the cloud cover. Additionally, the seasonality in mean EKE is comparable between observations (altimetry) and the model (not shown), which also suggests the limitations of the MUR dataset.
Other tracer-based front detection algorithms, such as the objective Cayulla Cornillon Algorithm (Cayula and Cornillon, 1995) and the Nieto method (Nieto et al., 2012), could have been used to detect temperature fronts in the MC, but they do not provide gradient magnitudes. Moreover, we did find a temperature front intensity threshold and spatial variability similar to what is documented in Nieblas et al. (2014), despite using a different front detection algorithm.
Finally, it would be interesting to use Lyapunov Exponents as in Hernańdez-Carrasco et al. (2012); Bettencourt et al. (2012) and Siegelman et al. (2020b;2020a) to bring a dynamic-based approach to our multidecadal tracer-based front database.

Spatial variability
We have leveraged the 3D component of the well-validated MOZ36 to go beyond the surface and explore temperature front activity in the MC ocean interior as well. To our knowledge, this is the first temperature front spatial and seasonal variability analysis extending to the ocean interior in the MC.
The mean maps of frontal intensity and probability (Figure 2) from MOZ36 and MUR SST present similar spatial patterns to what is shown in Nieblas et al. (2014). Additionally, the vertical distribution of temperature front intensity is consistent with the stratification profiles from MOZ36 and CARS2009 ( Figure S4). Indeed, the Islands (resp. Agulhas) region has the deepest (resp. shallowest) thermocline and frontal intensity peak, at about -200 m (-80 m). It is worth noting that as MOZ36 does not consider the run-off from the Zambeze river, we expect thermal shelfs fronts over the SofalaBank to be underestimated. In all 4 regions, front intensity is maximum at around -100 m, in the vicinity of the mean thermocline. The maximum intensity fronts at around -100 m may be explained by strain-induced vertical velocities at the edge of eddies, which intensify below the mixed layer and up to -400 m below it (Keppler et al., 2018;Siegelman et al., 2020a).
ADF develop at the border of large anticyclonic eddies, or Mozambique Rings, on the western side of the MC, as observed by Biastoch and Krauss (1999) and Schouten et al. (2003). On the other side of the channel, CDF form at the edge of smaller cyclonic eddies, mostly on the southeastern side of the channel and driven by the SMEC, as seen in Gründlingh's work (Gründlingh, 1995). In MOZ36, the temperature front signature of the anticyclonic Mozambique Rings seems to be limited to the upper layers of the water column. However, some discrete observations (de Ruijter et al., 2002) and MOZ36 modelled current (Miramontes et al., 2019) have shown that the dynamics of the Mozambique Rings could indeed reach the bottom. It is still unclear why the temperature front signature of these anticyclonic rings remains weak below 500 m while their associated currents have been previously reported at much deeper levels in the model.  Seasonal variability in the MC ceases to be significant below -150 m for the Islands subregion, below -350 m and -450 m for the NMC and SMC subregions (respectively) and below -70 m in the Agulhas subregion (Figure 3). This is aligned with the types of fronts that characterize each subregion, as summarized in Table 1. The NMC and SMC subregions are eddy-rich subregions generating mostly ADF and CDF fronts which are more likely to reach deeper in the water column than topographic fronts, and thus propagate surface seasonal signals at depth. We suspect this explains why the seasonal variability in the NMC and SMC subregions reaches deeper layers in the water column. Potential drivers of the seasonal variability for each region are summarized in Table 2.

Islands and NMC subregions
In the Islands and NMC regions, SST front intensity climatology shows two local maxima: a moderate one in May and a greater one in November (Figure 7). A similar semi-annual cycle at the Seychelles Chagos thermocline ridge (Hermes and Reason, 2008) suggests that the Seychelle-Chagos upwelled waters might be one of the possible drivers of temperature front seasonal variability at the surface of both northernmost regions. Additionally, the presence of "cool water" events near Angoche during August to March (Malauene et al., 2014), and upwelling along the Mozambican coast at 11-16°S during the North East monsoon (November-March) (Saetre and Da Silva, 1984) could explain the higher intensity of the second peak in austral spring and summer, especially in the NMC. Another driver behind the local winter maximum could be the elevated EKE in this area responding to the seasonal increase of the NEMC (Figure 9) (Obura et al., 2019).
In the Islands subregion, the -100 m front intensity is greater than at the surface and peaks in late austral autumn. The signature 2-peak of the surface signal has disappeared entirely, leaving only one clear maximum in May-June (Figure 7). This suggests that the -100 m fronts may not be driven by the same processes than the surface ones. Since most fronts in the Islands region belong to the topographic class, and since the NEMC peaks in austral winter (Obura et al., 2019), -100 m temperature fronts in the Islands region seem to be driven solely by seasonal variations of the NEMC and its subsequent interactions with the local topography.
In the NMC area, the -100m front intensity is higher than at the surface and peaks in austral summer, with a second local maximum in austral winter (Figure 7). Those two peaks can be explained by the types of front predominantly found in this region: topographic fronts, and especially ADF (Table 1). On the first hand, as seen in the Islands subregion, topographic fronts follow the variability of the NEMC, which is stronger in austral winter (Obura et al., 2019), explaining the winter peak of frontal activity. On the other hand, ADF are shaped at the edge of anticyclonic eddies forming in this region (Schouten et al., 2003;Halo et al., 2014) and are thus likely linked to variations of EKE, which would explain the frontal activity maximum in austral summer, coinciding with peaks of EKE (Figure 9) At -500 m and -1000 m, there is not much seasonal variation anymore and frontal activity weakens with depth ( Figure 9). This is consistent with our preliminary front typology since ADF and topographic fronts, which represent the majority of fronts in Islands and NMC regions, have been found to present a weak deep signature.

SMC and Agulhas subregions
The FFT periodograms of the two southernmost subregions show that the yearly signal is stronger than the biannual one ( Figure 6). The weaker biannual signal can be explained by the lesser influence of the monsoon below 20°S (Schott et al., 2009). The climatology of surface temperature front intensity in the  (Figure 7) shows stronger frontal activity during austral winter and spring which coincides with upwelling cells at the southern tip of Madagascar driven by the SEMC (Ho et al., 2004;Ramanantsoa et al., 2018) and the local winds (Lutjeharms and Machu, 2000;DiMarco et al., 2000). Frontal activity in the SMC and Agulhas regions also show a slight increase around December (Figure 7). One possible explanation is the Angoche upwelling cell, which brings cooler water to the surface in austral summer (Malauene et al., 2014). Another possible driver could be the eddy activity, which peaks in austral summer in these regions. Indeed, since the SMC and Agulhas regions seems to host mostly CDF and topographic fronts, it is likely that the temperature front variability in these regions depends on EKE variability.
In the SMC at -100m depth, the first local maximum in austral summer (Figure 7) coincides with a -100 m EKE peak (Figure 9), which suggests that frontal activity there is driven by MC eddies. The second peak in austral winter can be attributed to the local upwelling cells south of Madagascar island (Ramanantsoa et al., 2018), which was also the main driver of surface frontal activity in this region. In the Agulhas region, the -100 m front climatology (Figure 7) shows much weaker seasonal variation, also seen on the FFT analysis ( Figure 6), with a slight maximum in austral winter, when the South Madagascar upwelling cells are most active (Ramanantsoa et al., 2018). The absence of austral summer peak, like it is observed in NMC and SMC, can be explained by the comparably weaker EKE seasonal variations in the Agulhas region ( Figure 9). It is worth noting that the FFT analysis ( Figure 6, and Figure S6) does not show a clear annual or bi-annual period for the SMC and Agulhas regions from -100 m depth and below, which could be linked to the higher inter-annual variability observed in this region (Palastanga et al., 2006).
-500 m and -1000 m fronts in the SMC and Agulhas subregions are more intense than in the Islands and NMC at the same depth, possibly due to the specific typology of fronts most often encountered in the south MC (Table 1). In the SMC and Agulhas, the majority of fronts belong to the CDF category, which have a stronger depth signature, hence the greater -500 and -1000 m front intensity compared to regions where ADF dominate.
Finally, a modest peak at about 4-month can be identified in the FFT analysis periodograms (Figure 6), especially at the surface of the Islands, SMC and Agulhas subregions. This could be attributed to Indian Ocean basin resonance characterized by 120-day period, as detected by both idealized and realistic models (Han et al., 2011).

Main drivers of frontogenesis in the MC
We have thus documented thermal front seasonal variability in the MC. The next section addresses the relative importance of the main drivers of frontogenesis in the region. A clear distinction between those physical processes would require carrying dynamical diagnostics such as eddy kinetic budgets [e.g., Conejero et al. (2020)], eddy kinetic energy decompositions [e.g., Martínez-Moreno et al. (2019; or frontogenesis function decompositions [e.g., Koseki et al. (2019)] based on the model outputs, which would fall beyond the scope of this work. We nonetheless argue that thermal fronts in the MC depends on mean and seasonal circulation, seasonal upwellings and variations of EKE. To do so, we first define our seasonal division: In the Mozambique Channel monsoon reversals occur first in April to June and then in October to November, marking the separation between winter and summer (Collins, 2013;Obura et al., 2019). We thus define winter as the period from April to September, and summer as the period from October to March.

Mean circulation
We can infer the role of the mean circulation on thermal front intensity (especially on the quasi-permanent fronts since no seasonal variability is involved here) from mean maps of surface current (northward and eastward velocities) and SST generated by the MOZ36 configuration ( Figure S7). From the mean SST map ( Figure S7.B), we can observe a strong thermal front at the northern tip of Madagascar, following the path of the North East Madagascar Current (NEMC), also visible in the mean current map ( Figure S7.A) and the mean thermal front intensity (Figure 2). Similarly at the southwestern tip of Madagascar, a strong thermal front follows the path of the South East Madagascar Current (SEMC) toward the Agulhas current. Additionally, we can identify along-shelf thermal gradients of moderate intensity at the southeastern tip of Madagascar and off Maputo (Delagoa Bay), which are probably due to quasi-permanent coastal upwelling cells (Ramanantsoa et al., 2018), as well as a colder tongue extending along the outer parts of the Sofala Bank, which could be indicating the presence of a quasi-permanent shelfbreak upwelling there, as documented by Hill and Johnson (1974); Gibbs et al. (2000) and Rossi et al. (2010).

Seasonal upwellings
Seasonal upwelling cells such as the Seychelles-Chagos upwelling (Hermes and Reason, 2008), the Angoche upwelling (Malauene et al., 2014) and, to a lesser extent, the South Madagascar upwelling cells (Ramanantsoa et al., 2018) mentioned in this work are strongly localized in space and time. The Seychelles-Chagos upwelling turns on in austral winter (Hermes and Reason, 2008), while the Angoche upwelling is active during the August-March period (Malauene et al., 2014). The first upwelling core south of Madagascar is present all year round but reaches its maximum at the end of austral winter, while the second (much weaker) core is most active from October to January (Ramanantsoa et al., 2018). In the Islands and NMC regions, the mean thermal front intensity peaks first in May-June (Figure 7), which corresponds to the first active period of the Seychelles-Chagos upwelling system. This front intensity peak is mostly felt in the Islands subregion which is closest to the Seychelles-Chagos upwelling cell. We thus postulate that this first front intensity peak at surface and -100 m depth in both Islands and NMC subregions depends on the Seychelles-Chagos upwelling system. The second peak of frontal intensity in the same two northern subregions happens in November-December (Figure 7), and especially in the NMC region. We argue that this second peak is driven first by the Angoche upwelling and to a lesser extent by the Seychelles-Chagos upwelling. The seasonality of this second peak coincides with the Angoche most active period, which is located in the NMC, and the second (weaker) upwelling maximum in the Seychelles-Chagos located further upstream. The South Madagascar upwelling cells reach their maximum in winter, when the thermal front intensity peaks in both NMC and Agulhas subregions Figure 7). It thus suggests that the South Madagascar upwelling cells play an important role in the increased front intensity documented in the southernmost subregion during winter.

Seasonal current and EKE variability
Some quasi-permanent thermal fronts are found in the path of the NEMC, the SEMC and the Mozambique Current. The seasonal variability of these currents induces in turn variability of instabilities. We suggest that these instabilities, in turn, produce the sub to mesoscale activity that affects the intensity of those thermal fronts. The NEMC, and the EKE associated to its instabilities, are both stronger in winter ( Figure S8 and S9), which also corresponds to a higher thermal front intensity in the northern Islands subregion, especially at -100 m ( Figure S10 and S11). The Mozambique Current and the EKE in the western MC are stronger in summer ( Figure S8 and S9), which also corresponds to higher thermal front intensity in the western half of the MC, especially in NMC and SMC subregions ( Figure  S10 and S11).
Finally, the suspected main drivers of seasonal variability of thermal front intensity, for each season and each subregion, and at both surface and -100 m, are summarized in Tables 2 and  Table S1.

Conclusion and perspectives
We presented here the first thermal front seasonality analysis in the Mozambique Channel from the surface down to -1000 m depth. This work demonstrates that front dynamics in the MC are complex, and we suggest several possibilities to link temperature front seasonal variability to regional and remote ocean dynamics. Temperature fronts are generally more intense and more likely to occur over continental shelves, in the wake of the NEMC and SEMC, on the lee of islands, and at the edge of eddies. We show that temperature front intensity peaks below the mixed layer, between -100 and -200 m depth, which indicates that they are largely associated to aspects of the ocean dynamics rather than to mixed-layer processes. This consistent with the energetic mesoscale dynamics in this region.
Seasonal variations are felt down to -150 m in the Islands subregion, -350 m and -450 m in the NMC and SMC subregions (respectively) and -70 m in the Agulhas subregion. We argue that the main drivers of thermal front seasonal variability in the MC are upwelling cells (Angoche, South Madagascar, Seychelles-Chagos), current seasonal variability (NEMC, SEMC, and the Mozambique Current) and variations in EKE.
Since the MC circulation is connected to equatorial circulation in the Indian ocean, it is expected that temperature fronts experience fluctuations at interannual timescales. The Indian Ocean Dipole (IOD) seems to affect the variability of the transport (Ridderinkhof et al., 2010) and the SSH in the channel (Palastanga et al., 2006) through planetary wave propagation. A positive (negative) IOD event is linked to a weaker (stronger) NEMC (Ridderinkhof et al., 2010) which could potentially affect the intensity of fronts in the Islands and NMC subregions through fluctuations in current instabilities. Additionally, a change in SSH would also modulate the interactions between currents at the entrance or within the MC, which affects not only the EKE in the MC (Palastanga et al., 2006) but also the production of temperature front at the edge of eddies.
Another relevant path of investigation is suggested by Martıńez-Moreno et al. (2021) in their global multidecadal study on altimetry and SST satellite-based observations. Their analysis shows that, even though SST gradients in the northern part of the MC have generally weakened between 1993 and 2020, the EKE has significantly increased during the same period in the northern MC and decreased in the southwest MC. This multidecadal trend suggests that climate change could have a significant impact on eddy-driven frontogenesis in the channel, as temperature fronts could be weaker but more numerous. Li et al. (2020) have also demonstrated a significant increase of stratification over the past half-century in the MC, especially in the Islands and NMC subregions. This suggests a possible enhancing of frontal activity due to climate change, especially in the two northernmost subregions, yet to be studied.
Thermal fronts, within and below the mixed layer, are thought to create transitory biodiverse hotspots across all trophic levels. Our analyses of their variability sheds light on where and when they are the strongest, with each subregion behaving differently. Altogether, it provides critical insights about where, horizontally and vertically, it is more relevant to survey bio-physical interactions and to consider protection initiatives in the MC. Conservation stakeholders could indeed use our synthetic analyses of front spatial and seasonal variability to elaborate rules for an adaptive and dynamic management of both epipelagic and mesopelagic realms, as also suggested in Brito-Morales et al. (2020) work on climate velocities. More research is needed to better apprehend climate change implications for front variability and their associated biological impacts to design an effective and climate-resilient network of pelagic and coastal marine reserves in the MC.

Data availability statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions
FS with VR, BD, and VG designed the study. PP provided MOZ36 model output, and guidance for model analyses. JS retrieved the remotely sensed MUR SST data. FS produced the front database and subsequent analysis of the results. CM assisted with the High Performance Computing. VR, BD, and VG supervised the findings of this work and provided critical feedbacks. FS produced the figures and wrote the manuscript with input from all authors. All authors contributed to the article and approved the submitted version.