Internal tides in the central Mediterranean Sea: observational evidence and numerical studies

Internal tides are studied in the central Mediterranean Sea using observational data and numerical experiments. Both numerical results and observations indicate that the baroclinic variability in this area is dominated by the K1 diurnal tide. In agreement with previous studies, the diurnal internal tides have the characteristics of Kelvin-like bottom trapped waves. They are mainly generated by the interaction of the induced barotropic tidal flow with the steep bathymetric gradient connecting the Ionian Sea with the shallow Sicily Channel. The bathymetric gradient appears to be the major forcing shaping the propagation paths of the internal tides. The most energetic internal tides follow the steep bathymetric gradient, propagating southward and tending to dissipate rapidly. Other waves cross the continental shelf south of Malta and then split with one branch moving toward the southern coast of Sicily and the other moving toward the west. Internal tides propagate with a variable phase velocity of about 1 ms−1 and a wavelength of the order of 100 km. During their journey, the internal waves appear to be subject to local processes that can modify their characteristics. The induced vertical shear strongly dominates the vertical turbulence and generates vertical mixing that alters the properties of the water masses traversing the area. Barotropic and internal tides remove heat from the ocean surface, increasing atmospheric heating, and redistributing energy through increased lateral heat fluxes. Lateral heat fluxes are significantly greater in the presence of internal tides due to the simultaneous increase in volume fluxes and water temperatures.


Introduction
The Mediterranean Sea is a semi-enclosed anti-estuarine sea divided into two anti-estuarine sub-basins, the Western and Eastern Mediterranean, to the west and east of the Sicily Channel, respectively (Pinardi et al. 2015). Deep and intermediate water mass formation occurs in both sub-basins. The waters of Atlantic origin and the intermediate waters formed in the Levantine basin are present in both the eastern and western Mediterranean sub-basins. On their way through the Mediterranean, the newly formed Levantine intermediate waters and waters of Atlantic origin together cross straits and channels that shape the geometry of the Mediterranean Sea. Straits and channels are quite energetic areas in terms of currents intensities and play a fundamental role in the circulation of the Mediterranean Sea. In these areas, local processes modify the properties of the water masses (Béranger et al. 2004). The water masses emerging from these straits, channels or passages are fundamental in determining the Responsible Editor: Emma Young.
Plain language summary Although tides in the Mediterranean Sea are known to be small, in this study we found that their effects on the ocean interior can be crucial for driving ocean variability and modifying the characteristics of the water masses. Internal waves at tidal frequencies are generated in the central Mediterranean Sea by the interaction of the flow with the steep bathymetric gradient that connects the shallow area of the central Mediterranean Sea with the deep Ionian Sea. The internal waves generate vertical mixing that changes the properties of the water masses traversing the area, enhancing atmospheric heating, and redistributing energy through increased lateral heat fluxes.
1 3 preconditioning factors that control oceanographic processes in other areas of the basin. Investigating the dynamics in these passages clearly contributes to a better understanding of the circulation and variability of the Mediterranean basin.
In this study, we focus on the central Mediterranean Sea, specifically in the eastern Sicily Channel. This area is characterised by a complex bottom topography, with the shallow areas allowing only surface and intermediate waters to be exchanged. The circulation in this region has been extensively studied in the past (Lermusiaux and Robinson 2001;Astraldi et al. 1996Astraldi et al. , 1999Astraldi et al. , 2002Millot 1987;Lermusiaux 1999) and can be schematized as a two-layer exchange, with a surface flow of Atlantic Water from the west to east, and a deeper flow from east to west of water formed in the eastern basin mainly composed of Levantine Intermediate water. The surface circulation is characterised by jet streams forming meanders, vortices and eddies constrained and driven by topographic effects, internal baroclinic processes, and strong atmospheric forcing (Lermusiaux and Robinson 2001). The deep circulation is more stable, and it is mainly shaped by topography (Béranger et al. 2004).
The variability in this region, similar to the rest of the Mediterranean Sea, spans several temporal and spatial scales (Borzelli and Ligi 1999;Grancini and Michelato 1987;Manzella et al. 1990;Pierini 1996). Although tides are relatively small in the Mediterranean Sea, tidal effects can be locally significant at high frequencies or small scales (Grancini and Michelato 1987;Artale et al. 1989; Mosetti and Purga 1989;Molines 1991;Lermusiaux and Robinson 2001). In particular, while sea level tidal ranges are small and mostly dominated by the M2 semidiurnal component (Lozano and Candela 1995) tidal currents can be substantial and are dominated by the K1 diurnal component (Cosoli et al. 2015;Poulain et al. 2018). Askari (1998) and Artale et al. (1989) observed internal waves with diurnal periods in the Sicily Channel region and both interpret them as the result of the tidal forcing interacting with the bathymetric gradient (Baines 1982). Artale et al. (1989) found that the variance explained by the 24-h internal wave account for a large fraction of the total observed variance. The authors pointed out that diurnal frequencies at this latitude are sub-inertial and the observed diurnal internal waves must be topographically trapped internal Kelvin-like waves. The weakness of the semidiurnal internal waves is not due to the proximity of the M2 amphidromic points but rather, as noted by Molines (1991), to the fact that the forcing mechanism for internal waves comes from tidal velocities dominated by the diurnal component, rather than from the sea surface elevation.
The present study is motivated by a sea trial conducted by the NATO Centre for Maritime Research and Experimentation southeast of Malta in the fall of 2020. The objective of this trial, the Mediterranean Rapid Environmental Picture 2020 (MREP20), was to study small-scale oceanographic variability and to assess its impact on sound propagation (Poulain and Oddo 2020). The collected data have been recently analysed (Poulain et al. 2023) and the results confirm that the dynamics on the Malta Bank are dominated by a large amplitude tide whose main components are diurnal baroclinic oscillations and vertical fluctuations of the isopycnals below the seasonal thermocline. Their study confirms that tidal currents are more relevant than tidal elevation in forcing internal dynamics, the former being dominated by the diurnal component (K1) and the latter mainly by semi-diurnal oscillations (M2). Diurnal oscillatory signals are enhanced near the shelf break and have the characteristics of an internal shelf (or Kelvin-type) wave propagating with the shelf (coast) on the right. Poulain et al. (2023) estimated a phase velocity of about 1 ms −1 and a wavelength on the order of 100 km for the K1 internal wave. However, the observational results remain ambiguous, and the theory of internal waves can be applied with some difficulty to determine their properties and kinematics. For example, it is not clear where the observed diurnal internal tides are generated (Poulain et al. 2023).
To better elucidate the properties and dynamics of these internal waves, we examined the results of numerical simulations performed during the MREP20 observational period. While numerical studies of internal tides in the ocean are quite widespread (Holloway and Barnes 1998;Niwa and Hibiya 2001;Di Lorenzo et al. 2006;Tanaka et al. 2010;Osborne et al. 2011, to name a few), there are few numerical exercises in the literature that address the dynamics of the internal waves in the Mediterranean Sea. We analysed the output of a 3D primitive equation hydrostatic ocean model that realistically represents a broad spectrum of dynamics on the shelf and slope in the central Mediterranean Sea. The agreement between model results and observations demonstrates the adequacy of the modelling hypotheses and allows us to investigate the space-time characteristics of the internal tides and their impact on local dynamics. After the introduction, the experimental design used in MREP20 is illustrated both from observational and modelling perspectives. Section 3 presents and discusses the model results. Summary and conclusions are provided in Sect. 4.

Observations: the MREP20 experiment
The MREP20 sea trial was conducted in the central Mediterranean Sea southeast of Malta by the NATO Research Vessel Alliance between 27 October and 8 November 2020 (Fig. 1). The main objective of MREP20 was to study smallscale oceanographic variability and evaluate its impact on sound propagation. During the experiment several moorings 1 3 with oceanographic and acoustic instruments were deployed, along with gliders, Lagrangian drifters and profiling floats. Additional oceanographic data were collected using the instruments deployed by the ship. In this study, due to the lengths of the observational records, we used only a subset of the data collected by the moored oceanographic instruments to validate the model results. Table 1 provides an overview of the data used. A more detailed description of the available data and data analysis methods can be found in Poulain et al. (2023) and Poulain and Oddo (2020).

Ocean model
The ocean model implemented for the numerical experiments is based on Nucleus for European Modelling of the Ocean (NEMO) primitive equation model, version 4.04 (Madec and the NEMO team 2012). The model domain extends from 13.1° to 16.5° E in longitude and from 34.0° to 37.0° N in latitude (Fig. 1). Bathymetry is from the General Bathymetric Chart of the Oceans (GEBCO; Weatherall et al. 2015). The horizontal resolution of the model is approximately 1.125 and 1.5 km along the zonal and meridional directions, respectively. The vertical grid is discretized with 141 depth levels with partial steps near the bottom to ensure a proper representation of the bottom topographic features (Barnier et al. 2006). The bulk formulas from Large and Yeager (2004) are used to iteratively compute surface fluxes at each model time step, using atmospheric fields of air temperature, humidity, wind, cloud cover and atmospheric pressure derived from the Global Forecast System weather forecasting system operated by the National Weather Service   (Flather 1976) for barotropic velocities, a flow relaxation scheme for baroclinic velocities and a radiation scheme for active tracers (Oddo and Pinardi 2008). The flow relaxation scheme uses a dumping scale of 3 days, and the radiation scheme uses 3-and 10-day dumping scales in case of inflow or outflow, respectively.
The vertical turbulence coefficients were calculated using the Turbulent Kinetic Energy closure scheme (Gaspar et al. 1990;Madec et al. 1998). The model uses a splittime step technique to solve slow and fast ocean dynamics in a computationally efficient manner. The baroclinic and barotropic time-steps are 180 s and 3 s, respectively. Advection of active tracers (temperature and salinity) is resolved using a 4th order Flux Correct Transport scheme (FCT, Zalesak 1979) while a bi-laplacian operator along iso-neutral surfaces resolves horizontal diffusion. The momentum equation is solved numerically using a 3rd order Upstream Biased Scheme (UBS, Shchepetkin and McWilliams 2005) for the advective terms and a bi-laplacian operator along iso-level surfaces for the viscous terms. To maintain numerical stability and minimize numerical diffusivity, spatially variable horizontal diffusivity and viscosity coefficients were used. The values of horizontal diffusivity and viscosity are homogeneous within the model domain (0.8 × 10 6 m 4 s −1 and 1.7 × 10 6 m 4 s −1 , respectively), while they increase linearly near the coast and open lateral boundaries from a distance of 10 points of the model grid. These diffusivity and viscosity values in the model interior, in the proximity of the coast or lateral open boundary are based on the results of several sensitivity experiments performed during the model setup phase. The model is constrained at the surface by adding a feedback flux term based on the difference between the satellite level-4 sea surface temperature (SST) product from CMS (Buongiorno Nardelli et al. 2013) and the predicted temperature in the first model layer. The feedback coefficient is − 120 Wm −2 K −1 at night and gradually decreases to 0 Wm −2 K −1 at noon, to account for the use of SST measurements only at night in the level-4 SST product. Light extinction coefficients are calculated interactively using daily chlorophyll values observed by satellites (Lengaigne et al. 2007) available through the CMS portal.

Numerical experiments
To understand the nature and properties of the internal tides observed in the MREP20 dataset (Poulain and Oddo 2020), three different model simulations were performed. All the experiments were initialized and one-way nested with the CMS Mediterranean Sea operational model (Clementi et al. 2021) which has a horizontal resolution of about 4.5 and 3.5 km along the meridional and zonal directions, respectively. The experiments were initialized on 1 October 2020 and ran until 30 November 2020. Data at the lateral open boundary were taken as daily means from the CMS operational service, thus filtering out tidal fluctuations. A reference experiment was performed without introducing the tidal forcing at the lateral open boundary or activating the tidal potential term in the momentum equations of the primitive equation model. A second experiment was performed by superimposing the eight major constituents (Q1, O1, P1, K1, N2, M2, S2, and K2) of tidal elevations and velocities from TPXO (Padman and Erofeeva 2004) on the lateral open boundary conditions and adding the tidal potential to the surface pressure gradient calculations. To distinguish between the effects of diurnal and semidiurnal tides, a third experiment was performed using only the diurnal tidal components. Results of the third experiment are discussed in Sects. 3.4 and 3.5 while the remaining sections of the manuscript discuss the results from the simulation accounting for all tidal constituents.

Comparison with observations
The analysis of the results begins with the comparison between the observed and modelled power spectra calculated from the vertical displacement of the isopycnals at the mooring locations (Fig. 2). The model results were only considered to match the observational periods (Table 1). We selected a set of isopycnals that can be tracked over the entire time-series while optimizing the vertical resolution. The depth of selected isopycnal surfaces were determined and the time series of vertical isopycnal displacements were determined by subtracting the average depth of the respective isopycnal surface. The isopycnal vertical displacement was assigned to the mean depth of the isopycnals. Power spectra densities are calculated using Welch's method (1967) by dividing the time series into 3.5-day window segments with 1-day overlap and applying a Hamming window to each segment. The vertical coverage of the spectra is a consequence of the sensor depths and the methods used (isopycnal displacement does not work in the mixed layer). At the OM1 location, both the model simulations and the observations show that most of the energy is at low frequencies, especially in the deeper layers. In the observations and in the simulation accounting for the tidal forcing, there are internal baroclinic local maxima at the K1 and M2 frequencies, although the simulated vertical structure only partially resembles the observations. Without the tidal forcing, the model tends to produce relative maxima with a vertical baroclinic structure at the local inertial frequency (~ 20.6 h). Toward the shelf, at the OM2 site, the amount of energy contained in the tidal frequencies increases considerably. Unfortunately, mooring OM2 was only equipped with a few sensors between 135 and 300 m depth, so it is not possible to analyse the entire water column. However, similar to OM1, the model simulation forced with the tides shows significant maxima at both the K1 and M2 frequencies, while the simulation without tides shows maxima at the inertial frequency.
At the three shallow moorings (OM3, OM4 and WW), the presence of internal tides was evident at both K1 and M2 frequencies. However, the K1 frequencies are much more energetic than the M2 frequencies. The vertical extent of the maximum varies among the different sites due to the different depth and stratification of the water column. The simulation without tides does not reproduce the observed energy distribution and still shows the internal waves at the inertial frequency. With the introduction of the tides, the model captures the vertical energy structures observed at the different mooring locations quite well. Few notable differences exist for OM4 where most of the energy is trapped near the bottom in the model simulation and only one relative maximum occurs at the K1 frequency at mid-depth. Differences between model results and observations may be ascribed to several factors, most probably the misrepresentation of Fig. 2 Power spectra density computed from the isopycnal vertical oscillations at the different mooring locations. Y-axis is the mean depth of the isopycnals. Colours are the logarithm of power spectral density expressed in m 2 s −2 cph −1 . From top to bottom: OM1, OM2, OM3, OM4 and WW moorings. Columns a, b and c show observations, model results without tidal forcing, model results with tidal forcing, respectively bathymetry being the dominant one. A preliminary comparison (not shown) of bathymetric data collected during the sea-trial and the dataset used to create the model grid revealed a generally good agreement in the large-scale bathymetric features of the studied area but differences were observed in the exact isobath positions. In agreement with Poulain et al. (2023), we conclude that internal waves with diurnal frequency are enhanced on the shelf, while semidiurnal internal tides are also present at the base of the mixed layer but are generally less energetic than diurnal tides.
It is clear that by introducing the tidal forcing, the model is able to reproduce the dynamics of the internal tides, while without tides the local energy maxima occur at the inertial frequency. The vertical structure is only qualitatively similar between the model forced with tides and the observations. The differences may be due to errors in bathymetry, or in the vertical structure of the water column. However, it is clear that the internal tides and their spatial distribution are adequately represented in the model. In the rest of the paper, we will focus on the discussion of diurnal internal tides because they are the most energetic.

Temporal and spatial distribution
The temporal coverage of the observational records is relatively short (WW, the longest record is about 12 days), making it difficult to capture the temporal evolution of internal tides from the observations. However, given the good agreement between the modelled and observed energy distributions at the mooring locations, we can examine the time and space domains of the model to gain a better understanding of the local dynamics. Figure 3 shows the time series of the amplitude of the diurnal tidal current (extracted from the model simulation using the MATLAB® t_tide algorithm, Pawlowicz et al. 2002) along with the Morse wavelet coefficients (Olhede and Walden 2002;Lilly 2017) calculated from the isopycnal vertical displacement at the WW location and at the depth of maximum K1 energy (approximately 150 m in the observations and model simulation without tides, and 130 m in the model simulation accounting with tidal forcing, Fig. 2). Wavelet coefficients are shown for the two simulations and for the observations. In the model simulation that takes into account the tidal forcing, the temporal variability of the wavelet coefficient resembles the cycle of neap-spring tides. As expected, the energy at K1 or M2 frequencies is very weak in the simulation without tides. Although short, a temporal evolution of the wavelet coefficients can also be seen in the observational dataset, with the energy increasing during the study period in agreement with the occurrence of the spring tide. However neap-spring tides and oscillations in the amplitude of the internal tides are not perfectly in phase. This may indicate that internal tides are not generated exactly at the WW site. However, the model Following Nash et al. (2012), vertical velocity is used to identify internal tides, their distribution, and intensities. The spatial distribution of the vertical velocity energy associated with the internal K1 tides, in the entire model domain, is shown in Fig. 4 for three different 8-day time windows centered on the spring tide events. The vertical velocities are integrated in a 50-m thick bottom layer after removing the barotropic component (which accounts for barotropic tides, surface fluxes and bottom boundary conditions): where w ′ is obtained by first removing the barotropic component. Sensitivity experiments (not shown) indicate that the results are robust regardless the thickness of the bottom layer. The spatial distribution of K1 energy strongly resembles the bathymetric features, with higher energy associated with the largest bathymetric gradient (Fig. 5, left panel). The highest energies are near the bathymetric slope of the Ionian Sea and near the canyons west and northwest of Malta, with the lowest values observed on the Malta Bank. Temporal where u ′ and p ′ are the baroclinic velocity fluctuation and baroclinic pressure anomaly respectively, and angle bracket (< >) indicates the average over one, or more, tidal periods. Barotropic to baroclinic energy conversion rate (C topo ) is computed as: where u tb is the tide-induced barotropic velocity and H is the bathymetry of the model (Kurapov et al. 2003). Assuming that the local time derivative and the advection of diurnal internal tide energy terms are small, dissipation ( ) of internal tides is estimated as: To compute u ′ and p ′ in Eqs. 2 and 3, model baroclinic velocities and densities were bandpass filtered to retain only variability at K1 frequency (± 10% of K1 period, which is 23.93 h). In Fig. 5, the three terms (Eqs. 2, 3 and 4) averaged between 2 and 10 November are shown. In a, the tidal ellipses extracted from the model simulation at K1 frequency are overlaid with the topographic energy conversion term. Similar to the energy distribution, F , C topo , and are also characterized by very similar pattern and values during the other two events (not shown). Large positive C topo values are present on the Ionian bathymetric slope and south of the Malta Bank approximatively following the 300 m isobath (Fig. 5a). Negative C topo values, indicating energy converted from baroclinic to barotropic motion, south of Malta coincide with the relative minima observed in the energy distribution (Fig. 4b). The tidal ellipses show that in areas where the bathymetric features are less pronounced or where the major axes of the tidal ellipses are parallel to the slope, conversion is less effective. We conclude that topographic energy conversion is one of the most important mechanisms generating diurnal internal tides in the area and also determining their variability. The vertically integrated flux vector ( Fig. 5b and sketched in Fig. 1b) clearly indicates that the interaction of the barotropic flow with the bathymetry generates two different waves and wave patterns which are trapped by the topographic slope and propagate along this slope with the shelf to their right. The dissipation term distribution (Fig. 5c) indicates that south of Malta part of the energy removed from the internal tide is lost. Energy dissipation near the generation sites is mainly associated with the dissipation of internal tides having higher-order vertical modes .
The zonal and meridional Hovmöller diagrams of w it band-pass filtered at the K1 frequency are shown in Fig. 6. The band-pass filter used is a zero-phase 4th order Butterworth filter with a lower and higher cutoff periods of 21.537 h and 26.323 h, respectively (± 10% of K1 period, which is 23.93 h). Velocities are masked according to the energy distribution shown in Fig. 4b, to ensure that only regions with sufficient associated energy are considered. It is difficult to discern a clear pattern in the zonal direction. A slightly coherent pattern is observed from 15.5° to 16.0° E where the internal waves tend to propagate westward. West of 15.5° E the slope of the contours fluctuates, with the spatial scale becoming smaller toward the west. The meridional section is characterized by a more consistent pattern for the Fig. 6 Hovmöller diagrams of w it band-pass filtered (± 10%) at K1 frequency. a Meridional average (X axis are longitudes). b Zonal average (X axis are latitudes) most part. South of 36.0° N, the waves tend to propagate southward. North of 36.0° N the wave moves northward but the field appears quite discontinuous, with a small portion of the section where the direction of propagation fluctuates. Figure 7 shows maps of w it , band-pass filtered at K1 frequency (± 10%), during a particular event. Again, velocities are masked according to the energy distribution. Snapshots every 3 h are shown from 00:00 to 9:00 UTC on 7 November. The propagation path can be seen more clearly. Large w it values are associated with the slopes that connect the Ionian Sea with the Malta Bank, where topographic energy conversion is also large (Fig. 5a). The modelled wave pattern is consistent with the vertically integrated energy flux vector (Fig. 5b) and the interpretation of the Kelvin wave trapped near the bottom (Poulain et al. 2023;Artale et al. 1989). The internal wave moves while maintaining the bathymetric slope to its right. The most energetic wave packets move southward along the slope and clearly interact with the complex bathymetry. Other wave packets manage to pass the slope and move westward south of Malta. West of Malta, two propagation branches are evident. One branch moves northward toward Sicily, and a second branch continues its journey westward.

Single event analysis
We analyse the two identified main waves in the model results individually, by considering the sections above the isobaths of 320 m and 700 m (Fig. 1b, red and green lines,  respectively). The band-pass filtered integrated vertical velocities, w it , are sampled along the two sections and the corresponding Hovmöller diagrams are shown in Fig. 8. For each point along the section, the model results falling within a 4-km radius (2-3 model points) are averaged so that the velocity intensities are smoothed and are graphically shown normalized by their maximum value to aid interpretation.
Following the isobaths appears to be an efficient method for tracking the internal waves as the Hovmöller diagrams show more coherent patterns with respect to those shown in Fig. 6. Overall, the diagrams confirm the previous interpretation Fig. 7 Snapshots of w it band-pass filtered at K1 frequency (± 10%) masked only in the regions of high energy. Bathymetric contours are overlaid that the waves propagate from the continental slope connecting the deep Ionian Sea to the relatively shallow Sicily Channel. The waves propagate with variable phase speed westward in the red transect and southward in the green transect consistent with the bathymetric gradient.
Along the "red" transect ( Fig. 8, top panel), about 140 km from the origin, a clear discontinuity of velocities is observed (in agreement with the energy minimum, Fig. 4). This corresponds to the turning point west of Malta, where waves deviate and begin to move northward. Given the coherence between the previous and subsequent points along the section, we assume that the main signal is due to the propagation of the waves from the eastern sector of the domain. However, geometry and local energy transfer (Fig. 5 a and c) may play an important role here. Between 250 and 290 km (see Fig. 1 for a graphical localization), the phase propagation oscillates between positive and negative values. Here, the waves change from a southward to westward motion and the bathymetric gradient has also some discontinuities (Fig. 1b). We argue that following the isobath in this particular area is not an accurate method.
Along the "green" transect, we see that the energy is much lower in the southern sector, far from the point of origin (Fig. 8, bottom panel). In this region, C topo has large negative values, subtracting energy from the baroclinic component. Although a general behaviour is easily observed, the waves do not propagate with constant phase speed along the section. Between 100 and 130 km from the origin of the section, there is no clear slope in the Hovmöller. In this region, there is a meander in the 700 m isobath and dissipation occurs (Fig. 5c).
To quantify the wave properties, we calculated the cross power spectrum between successive points along each transect. Similar to the w it calculation, the cross-power spectra are calculated by dividing the time series from 4 to 12 November into 3.5-day window segments with an overlap of 1 day. From the phase shift, Δ , and from the distance between two consecutive points, d, we calculated the wavelength ( = 2 d∕Δ ) and the phase speed ( , where is the frequency of the diurnal internal wave). Given the adopted method, the estimated quantities refer to the relative orientation of the two considered points. The results are shown in Fig. 9.
The two sections differ in mean properties and scales of variability. The phase speed of the "red" section ( Fig. 9, top panel) oscillates between 0.5 and 1.6 ms −1 in absolute values, corresponding to a wavelength between 50 and 130 km. These values derived from the model agree well with the those estimated from the observations (Poulain et al. 2023). As the waves turn around Malta, they reach a minimum speed and wavelength of about 0.5 ms −1 and 40 km, respectively. Between 180 and 200 km from the origin, a change in the sign of the phase speed is observed, which is less evident in the Hovmöller diagrams. In this area, south of Malta, the bathymetry becomes more complex and exhibits a relatively deep trench parallel to the wave propagation direction and the bottom conversion term also has positive values. We speculate that the nonlocally generated wave may be undergoing local energy exchange. The oscillations observed in the Hovmöller diagram between 250 and 290 km correspond to changes in the sign of the phase speed. Fig. 8 Hovmöller diagrams of normalized w it band-pass filtered at K1 frequency (± 10%) along red and green transect shown in Fig. 1. The X axis follows the transect path (red transect in the top panel, green transect in the bottom panel). Cumulative distance from the transect origin is provided in the bottom X axis labels, while geographical coordinates are provided in the top X tick labels. Y axis indicates time in dd/mm format 1 3 The "green" section ( Fig. 9, bottom panel) is characterized by relatively small values of phase speed and wavelength in its southern sector (corresponding also to low energy values, Fig. 8). Approaching the northern part of the section, the phase speed and wavelength have values similar to those in the "red" section, with a speed of about 1 ms −1 and a wavelength of about 90 km. The change in the sign of the phase speed generally corresponds to the meandering of the 700 m isobath. We argue that we do not exactly follow the internal wave propagation path along the 700 m isobath, and that the values in certain locations are representative of the main propagation axes (long and fast waves), while in other areas, they are almost perpendicular to the propagation path (short and slow waves).

Impact on vertical mixing
The importance of mixing associated with the shear and breaking of internal waves in the general ocean circulation and energy budget has been extensively studied (Ferrari and Wunsch 2009;St. Laurent and Garrett 2002). The turbulent closure model used to compute the vertical eddy viscosity and diffusivity coefficients is based on a prognostic equation for the turbulent kinetic energy, and a closure assumption for the turbulent length scales. The time evolution of the turbulent kinetic energy results from the balance between a positive term related to the vertical shear, and negative terms related to stratification, vertical diffusion, and dissipation. Thus, in the present model configuration, there is no specific parameterization of mixing associated with internal wave breaking. However, the presence of internal waves feeds the shear production term in the prognostic equation for the turbulent kinetic energy term, resulting in higher vertical mixing coefficients. To understand and quantify the effect of the increased vertical mixing, we begin by comparing the vertical turbulence coefficients with and without tides along the two transects and averaging the model results from November 4 to 12. Figure 10 shows the logarithm of the ratio between the vertical mixing coefficients predicted in the simulation with and without tidal forcing. As expected, internal tides increase the vertical turbulence in both sections, especially near the bottom. Along the "red" section, the increase in vertical mixing is quite homogeneous from 200 m downward. It reaches peak values where the presence of internal tides increases the vertical turbulence coefficients by more than 3.5 times compared to the simulation without internal tides. There is only one region where the weekly averaged vertical turbulence coefficient is lower in the presence of internal tides, between 160 and 220 km from the origin of the section. This also corresponds to the region where the wave propagation direction has the opposite direction (Fig. 9), confirming the specificity of the area and the fact that the method of following the isobath is likely to fail.
Similarly, the vertical mixing coefficients along the "green" transect in the simulation with tidal forcing are significantly larger than those without the tides. There are two areas where the logarithmic values are higher than 1, between 60 and 90 km and between 120 and 140 km from the origin of the transect. Both areas are characterized by the highest energy content. The first section corresponds to an area where the internal tides are well traced following Fig. 9 Phase speed (Y axis) and wavelength (color-coded) of the K1 internal tide along the red (top) and green (bottom) transect shown in Fig. 1 the isobath, the phase speed and the wavelength are quite homogeneous (Fig. 9) and the maximum associated energy is observed (Fig. 8). The second section shows high variability in terms of phase speed and wavelength.
In both sections, there are also areas where the vertical mixing at the bottom of the mixed layer is noticeably increased. To distinguish between the effects of semi-diurnal and diurnal internal tides in the local dynamics, the vertical Fig. 10 Logarithm of the ratio of weekly (4-12 Nov) averaged vertical turbulent viscosity coefficients along transects in Fig. 1 (red top, green bottom), with and without tides. Positive values indicate that the simulation with internal tides has larger vertical mixing Fig. 11 Maps of the logarithm of the ratio of weekly (4-12 Nov) averaged vertical turbulent viscosity coefficients. a Integral between surface and 75 m depth between the simulation accounting for all tidal constituent and simulation without tides. b Integral between 75 m depth and the bottom between simulation accounting for all tidal constituents and simulation without tides. c As panel a but between simulation with all tidal constituents and simulation with only diurnal tides. d As b but between simulation with all tidal constituents and simulation with only diurnal tides turbulence coefficients of the three model simulations are compared in Fig. 11. For the comparison, the water column was divided into two layers: the first layer is obtained by integrating from the surface to 75 m depth and is representative of the mixed layer; the second layer, integrated from 75 m depth towards the bottom is representative of the bottom layer. The logarithm of the ratio of the vertical turbulent coefficients indicates that the semi-diurnal component plays a minor role in the vertical turbulence in both layers. In the mixed layer, the introduction of the tides results in an overall reduction in vertical mixing (Fig. 11a). This is due to the increased stratification as a result of the increased surface heat gain deriving from the introduction of tides ( Fig. 12a and associated discussion in Sect. 3.5). Only in the areas surrounding Malta, the near-surface vertical mixing is reduced. This is probably a consequence of surface current shifts due to the effect of the tides. In the large area north of Malta, the change in vertical mixing is probably due to the westward shifting induced by tides of the Atlantic Ionian Stream (Lermusiaux and Robinson 2001) and associated water masses. Although vertical mixing in the surface layer is almost halved in large areas of the model domain by the introduction of tides, it remains sufficiently large to ensure layer homogeneity. Below the mixed layer depth, the increased vertical mixing deriving from the introduction of the tides strongly follows the bathymetric features. In the deeper areas of the model domain, vertical mixing is unaltered by the introduction of the tides. The semi-diurnal tidal component clearly plays a minor role in vertical mixing. Both in the surface and deep layers, positive and negative values of the logarithmic of the vertical turbulence coefficient ratio alternates geographically. Spatial scales in the surface layer, suggest the effect of semi-diurnal tides is to shift the large-scale currents, while in the deeper layers the semi-diurnal tide effect is more on sub-mesoscale features distribution.

Heat balance and impact on the water mass characteristics
The overall effect of barotropic and baroclinic tides throughout the model domain is summarized in Fig. 12. The enhanced vertical mixing acts in the direction of removing heat from the surface to warm the deep layers. Despite the relaxation to the satellite SST, the reduced surface temperature and the larger surface currents generate an increase in downward heat flux, which is due to the bulk formulation we use to calculate the air-sea fluxes. The top panel in Fig. 12 shows that the total surface heat flux difference between the simulations with and without tides is constantly about 2 Wm −2 after 10 days of spin-up. The surface heat flux difference between the simulation with all tidal components switched on and that with only the diurnal components is negligible. The surface heat flux is calculated as: where A is the area of our model domain, Q s , Q l , Q h , and Q e are the incoming shortwave radiation, the net longwave radiation, the sensible and the latent heat fluxes respectively, all computed according to Large and Yeager (2004). This means that the ocean is absorbing heat from the atmosphere to warm its interior. In their barotropic and baroclinic forms, tides do not significantly alter the net volume transport across the open lateral boundaries (Fig. 12b). Their main effect is to introduce an oscillation of the volume flux at semidiurnal tidal frequency. However, the flows across each of the open lateral boundaries are equally amplified by the tides and compensate for each other. This is especially true for the eastern and western open lateral boundaries corresponding to the main tidal propagation direction. Recalling that the mean transport in the intermediate and deep layers is from east to west (Lermusiaux and Robinson 2001), to understand the effects of the introduced dynamics on the water mass properties we analyse the tidal effect on the outflow at the western open boundary of the domain. The last three bottom panels of Fig. 12 show, respectively, the western boundary outgoing heat flux, the averaged temperature of the outgoing water, and the outgoing volume flux computed as follows: where all the integrals and A, the area of the transect, are computed at each time-step only where the zonal velocity is negative (outgoing). The addition of tidal forcing causes an increase in lateral heat flux at the western boundary ( HF w ) , which varies at diurnal frequency. The increase in lateral heat flux is the combination of two factors: the increase in horizontal volume flux due to barotropic and baroclinic tides ( V ow ) , and the increase in the deeper layers temperature ( T ow ) due to the increase in surface fluxes and vertical mixing.
Analysing the horizontal maps of baroclinic transport integrated from 150 m depth to the seafloor and averaged over the period 4-12 November (Fig. 13a), we note that the introduction of the tides does not change the weekly averaged circulation characteristics in the area (Fig. 13b). In addition to differences in the positions of small structures mainly in the interior of the model domain, the main differences between the two simulations are in the narrowing of the strong inflow along the eastern boundary and the intensification of the outflow in the north-western sector of the model domain (Fig. 13b). On the other hand, the effects of internal tides on horizontal transport and vertical mixing become clear when comparing the vertically integrated temperature and salinity (from 150 m depth to the bottom). The bottom panels of Fig. 12 show the differences in temperature and salinity averaged over the same temporal window and integrated from 150 m depth to the sea floor, between simulations with and without tidal forcing. It can be seen that the north-western area of the domain is characterized by warmer temperatures, while the south-eastern region has lower temperatures and higher salinity. The first area is downstream of the main propagation path of the internal wave. It corresponds to the main outflow region and includes waters that were modified by the induced internal tidal mixing as they crossed the model domain from east to west. The changes in the eastern sector of the modal domain are probably due to the interaction with deeper water.
Finally, to examine the vertical structure of the differences, we chose a meridional section within our model domain (Fig. 1) and compared temperature and salinity characteristics averaged over the same week. The differences in temperature and salinity along the section are shown in Fig. 14. Model results are presented below 150 m, for clarity.
Along the section, after the water from the Ionian Sea has been transported through our domain and modified by the internal waves and their relative mixing, the differences seem to have a clear structure. The most striking feature is the minimum/maximum temperature difference between 200 and 400 m between 35.5° and 36.0° N. In this area, the waters circulate due to the presence of two interconnected anticyclonic eddies (Fig. 13a), which are slightly modified by the tides. North of 36°, there is a large area with much warmer temperatures and lower salinities between 150 and 300 m depth. This area corresponds to the internal waves that we have drawn with the red section. Thus, the internal waves play a significant role in changing the properties of the water masses. The correlation between temperature and salinity anomaly varies along the section and with depth according to the different water masses characteristics present in the area.
In the upper layers, surface waters locally modified (warm and salty) and waters of Atlantic origin (generally colder and fresher) interact and form complex horizontal and vertical structures, with Atlantic waters usually confided in a thin layer at about 150 m depth. Warm and salty waters of Mediterranean origin occupy the deeper layers. Fig. 13 a Baroclinic transport streamlines and intensity (colour) integrated from 150 m depth to the sea floor, for simulations with tidal forcing. b Streamlines of differences between baroclinic transport with and without tides with superimposed the differences in their modules. Temperature (c) and salinity (d) differences integrated from 150 m to the sea floor between model simulations with and without tides. All results are averages between 4 and 12 November

Discussion and conclusions
Inspired by recent analysis of data collected during a sea trial conducted by the NATO Centre for Maritime Research and Experimentation southeast of Malta in the fall of 2020 (Poulain et al. 2023), we conducted numerical simulations in the same area and same time period with particular emphasis on internal tide dynamics. Observational studies have revealed the presence of large amplitude tides whose main component is a diurnal baroclinic oscillation. However, some aspects of the internal tide dynamics, such as the generation site or propagation path, remain unclear. Surprisingly, despite the Mediterranean Sea circulation and dynamics have been largely investigated in the past, the relevance of internal waves (including tides) is poorly covered in the present literature. Considering that the studied area connects the Mediterranean sub-basins, a better knowledge of its dynamic clearly contributes to a better understanding of the overall Mediterranean circulation and its variability. Thus, in the present study, we characterized the internal tides and their effects on local dynamics by comparing numerical simulations performed with and without tidal forcing, or with only diurnal tidal components. In particular, we investigated the effects of the baroclinic tide on the surface and intermediate Mediterranean waters mixing.
Realistic and physically consistent simulations performed with a state-of-the-art high-resolution (about 1.3 km in the horizontal) model based on NEMO (Madec and the NEMO team 2012) provided new insights into the local dynamics. A reference numerical experiment was performed without introducing the tidal forcing and filtering out the tidal frequencies from the external lateral open boundary data; a second numerical experiment was performed by superimposing the eight major constituents of tidal elevations and velocities on the lateral open boundary conditions and also adding the tidal potential to the surface pressure gradient calculation; finally an experiment with only diurnal tidal component has been performed. Although the effect of diurnal and semidiurnal components is not purely additive, the latest experiment helps identify qualitatively the origin of variations in surface or bottom waters.
The ability of the model to reproduce the observed internal tide dynamics is evaluated in the observational space by qualitatively comparing the model results in terms of power spectra of isopycnal vertical displacement (Fig. 2). Although the simulated vertical structure only partially resembles the observation, it is clear that the internal tides and their spatial distribution are adequately represented in the tide-forced model simulation. Internal waves at diurnal frequency are enhanced on the shelf, semidiurnal internal tides are also present at the base of the mixed layer but are significantly less energetic than diurnal tides. Higher amplitudes of the diurnal baroclinic tides with respect to the semi-diurnal component agree with previous studies in the Mediterranean Sea on barotropic tides (Lozano and Candela 1995;Cosoli et al. 2015;Poulain et al. 2018).
At the WW site, the temporal variability of the wavelet coefficient associated with the maxima of the internal tides, resembles the neap-spring tidal cycles (Fig. 3a). The small phase shift between the energy maxima of the internal tides and neap-spring tides could indicate that the internal tides are not locally generated.
The spatial distribution of energy associated with the internal K1 tides, derived from vertical velocities, strongly resembles the bathymetric features, with higher energy associated with the largest bathymetric gradient (Fig. 1b) near the slope of the Ionian Sea, where large values of the topographic energy conversion term (Kurapov et al. 2003) are also observed. Large values of topographic energy conversion are also simulated in the southern boundary of the Malta Bank where barotropic diurnal tidal currents and bathymetric slope are orthogonal. In general, topographic energy conversion appears to be the most important mechanism for the generation of diurnal internal tides. Energy dissipation values near the formation site, usually associated with high-order vertical modes , agree with recent findings (Poulain et al. 2023) that the generated trapped diurnal internal waves can be described as a sum of several modes. Under these circumstances, a description of the propagation and kinematics of the waves in complete detail is difficult to achieve (Nagai and Hibiya 2022) and is beyond the scope of the present manuscript. However, the analysis of the internal tides characteristic provides robust results suggesting that only few low-order vertical modes propagate in the region.
Meridional and zonal Hovmöller diagrams of vertical velocity (Fig. 6) illustrate the complexity of the internal wave paths. Analysis of vertical velocity maps during a single event largely supports the interpretation of the model results. Consistent with the bottom-trapped Kelvin wave theory (Artale et al. 1989), the internal wave generated in the model at the Ionian Sea slope moves while maintaining the bathymetric slope to its right. While the main branch moves southward along the slope, part of the energy manages to pass the slope and move westward south of Malta where meets other locally generated internal tides. Part of the energy move southward (Fig. 5b), generating another local maximum approximately at 35° N 15° 30′ E (Figs. 4 and 5c). West of Malta, the internal wave further splits into two propagating branches, the most intense moving northward toward Sicily, and one moving westward.
As they travel, the waves propagate with variable phase speed and experience local intensification or attenuation as a result of bottom, mesoscale and mean current interactions. The phase speed oscillates between 0.5 and 1.6 ms −1 , corresponding to a wavelength between 50 and 130 km. Estimated values generally agree with the recent findings calculated directly from the MREP20 observations (Poulain et al. 2023). Below the mixed layer depth, the presence of internal tides also affects the vertical mixing, whose coefficients in large areas are dominated by the vertical shear term in the turbulence energy prognostic equations. In the steep area of the Ionian slope, vertical turbulence can increase in the model results by an order of magnitude as a result of internal tide dynamics. Differently from previous studies (Musgrave et al. 2016) semi-diurnal tidal components and induced baroclinic processes play a minor role in the local dynamics and in enhancing the vertical mixing. The introduction of the tidal forcing has the general effect of increasing the surface heat fluxes as well (despite the SST relaxation toward satellite observations); indeed, the simulation including the tidal dynamics receives more heat from the atmosphere. This increased surface flux in combination with the increased vertical mixing warms the deep layers modifying the water mass characteristics. In addition, barotropic and baroclinic tides also enhance east-west horizontal transport and increase exchange with surrounding areas. As a result, lateral heat fluxes are greatly increased in the simulation that accounts for tidal dynamics, due to both increased volume transport and to increased water temperature. The model results also suggest that internal tides and induced vertical mixing are a major driving factor in modifying the properties of the water masses traversing the area. Unfortunately, the MREP20 and other available climatological data do not have sufficient spatial and temporal resolution to support our findings on the effect of internal tides on the average proprieties of locally modified water masses. On the other hand, the breaking and dissipation of internal tides is not explicitly parametrised in our model, so the model results probably underestimate the importance of these processes for local dynamics.
The main mechanisms of internal tide generation and their impact on local dynamics have been investigated using realistic model simulations. Uncertainties remain, and a full description of the proprieties of internal tides will probably require more idealized model studies and further observational data.