Superstatistical analysis of sea surface currents in the Gulf of Trieste, measured by high-frequency radar, and its relation to wind regimes using the maximum-entropy principle

Two years (2021–2022) of high-frequency-radar (HFR) sea surface current data in the Gulf of Trieste (northern Adriatic Sea) are analysed. Two different timescales are extracted using a superstatistical formalism: a relaxation time and a larger timescale over which the system is Gaussian. We propose obtaining an ocean current probability density function (PDF) combining (i) a Gaussian PDF for the fast fluctuations and (ii) a convolution of exponential PDFs for the slowly evolving variance of the Gaussian function rather than for the thermodynamic $\mathit{\beta }=\mathrm{1}/{\mathit{\sigma }}^{\mathrm{2}}$ in a system with a few degrees of freedom, as the latter has divergent moments. The Gaussian PDF reflects the entropy maximization for real-valued variables with a given variance. On the other hand, if a positive variable, as a variance, has a specified mean, the maximum-entropy solution is an exponential PDF. In our case the system has 2 degrees of freedom, and therefore the PDF of the variance is the convolution of two exponentials.

In the Gulf of Trieste there are three distinct main wind forcing regimes: bora, sirocco, and low wind, leading to a succession of different sea current dynamics on different timescales. The universality class PDF successfully fits the observed data over the 2 observation years and also for each wind regime separately with a different variance of the variance PDF, which is the only free parameter in all the fits.

1 Introduction

Earth-observing systems provide us with an ever-increasing amount of data. This allows us today to move from considering low-order moments of fluctuating observations, averages and variances, to their probability density functions (PDFs), which give a complete characterization of the statistics and help expose the underlying dynamics. These PDFs are mostly non-Gaussian with fat tails, showing that extreme events are more frequent as compared to Gaussian statistics. Furthermore, when the system is subject to slow external forcing, the PDF of the fast fluctuations, characterized by the signal decorrelation time, evolves on a slower timescale. Examples of such dynamics are fast weather statistics (days) under slowly evolving climate conditions (tens of years) or fast sea surface current fluctuations (hourly) during slowly varying synoptic weather conditions (a few days), which is the problem considered here. Superstatistics, introduced by Beck and Cohen (2003), Beck (2004), and Beck et al. (2005b), is a formalism for considering dynamics on two separated timescales, as it considers a fast Gaussian process whose variance evolves on a slower timescale. This leads to long-term PDFs which are typically fat-tailed.

Figure 1 (a) Gulf of Trieste location (red rectangle) in the Adriatic Sea. (b) Gulf of Trieste zoom with the percentage of available high-frequency-radar (HFR) data (in multiple colours) in the selected period and principal axis ellipses (in red) of the $\mathit{\delta }=\mathrm{8}\times \mathrm{30}\min =\mathrm{4}\mathrm{h}$ velocity increments. The HFR baseline is shown with the red line between Aurisina and Piran. The HFR “P” grid point is shown with the black asterisk, and the closest WRF grid point is marked with the blue star and is called “WRF_P”. The selected HFR P grid point shows a high data coverage and is not close to the HFR baseline.

Superstatistics is a method widely used today in various scientific fields: environmental science (Weber et al., 2019; Schäfer et al., 2021), biology (Costa et al., 2022), statistical mechanics (dos Santos et al., 2022), and quantum science (Okorie et al., 2023). Superstatistics is applied here to the Gulf of Trieste, which is used as a natural laboratory to explore the air–sea interaction. The Gulf of Trieste is a shallow semi-enclosed basin in the northern Adriatic Sea (Mediterranean Sea; Fig. 1a): its maximum depth is approximately 25 m and its shape is a roughly 20 km × 20 km square opened on the western side. The Gulf of Trieste sea surface base circulation is driven by the basin-wide Adriatic cyclonic (anti-clockwise) circulation and by the internal thermohaline properties of the basin (Cosoli et al., 2012; Querin et al., 2021), showing a mean surface water outflow (Bogunović and Malačič, 2008). However, this base state is highly varied by two different forcings: the northern freshwater input from the Isonzo (Soča) River and the wind forcing (Querin et al., 2006; Cosoli et al., 2012; Cosoli et al., 2013; Querin et al., 2021). Moreover, the Gulf of Trieste and the northernmost part of the Adriatic shelf are also sites of North Adriatic Dense Water formation due to shelf convection (Jeffries and Lee, 2007; Pullen et al., 2007).

The main wind regimes consist of bora and sirocco wind events. The bora is a east-north-easterly katabatic wind with gusts reaching up to 50 m s⁻¹ and bringing cold and dry continental air from the north-east over the sea (Poulain and Raicich, 2001). The bora intensity is highly affected by the land orography and blows with particular strength over the Gulf of Trieste. Bora winds are most common during the cold season, but they can also occur during the summer. When the bora blows, it forces the surface water out of the gulf, and a replenishing bottom counter-current entering the Gulf causes an upwelling on its closed side (Querin et al., 2006; Malačič et al., 2001). Furthermore, bora events are very efficient in causing simultaneous cooling and evaporation by bringing cold and dry continental air in contact with the sea surface (Poulain and Raicich, 2001; Raicich et al., 2013). In contrast, the sirocco is a southerly moist and warmer wind channelled by the Adriatic coastal mountains that tends to occur year-round without a favoured month or season (Poulain and Raicich, 2001). The relatively smooth sirocco wind causes a sea-level rise in the northern Adriatic which drives an intense southward return flow when wind forcing relaxes (Cosoli et al., 2012).

High-frequency radar (HFR) is a powerful technology to measure horizontal surface currents and sea waves over a grid in a wide area of the sea. The HFR technology is founded on the principle of Bragg scattering of the electromagnetic radiation over the rough conductive sea surface (Crombie, 1955). Currents are obtained from the Doppler shift of radio waves backscattered by surface gravity waves of half of their electromagnetic wavelength. Each radar site allows us to estimate the component of the current towards or away from the receiving antenna (radial current component), and thus two sites are needed to obtain total currents. The distance to the backscattered signal is determined by range-gating the returns. Depending on the methodology used to determine the incoming direction of the scattered signal (also called “bearing determination”), commercial HFR systems can be divided into two main types: beam forming and direction finding. A detailed description of the radar technology and its uses can be found in Lorente et al. (2022).

This study aims to investigate the ocean currents in the Gulf of Trieste, focusing on the fluctuations of the recorded signals. This is done by adopting a superstatistical approach (Beck et al., 2005b). Superstatistics in fact is a method based on the development of statistics of statistics, which is well suited to the case of the Gulf of Trieste, where different sea current dynamics on different timescales superpose. Many studies (e.g. Ghil et al., 2011, and references therein) use the maximum-entropy principle as a tool in environmental science. Here, according to the Jaynes (1957) point of view, the entropy maximization will be applied to find the least-biased PDFs.

The paper is organized as follows: in Sect. 2 the used data are presented, in Sect. 3 the adopted methodology is explained, in Sect. 4 the results are shown, and finally the conclusions are in Sect. 5.

2 Data

Two classes of data sets are analysed: the sea surface horizontal HFR current data and the modelled wind data within the area of the Gulf of Trieste (northern Adriatic Sea, Fig. 1). The available data sets cover a time range of almost 2 years from 1 January 2021 to 18 October 2022.

2.1 The HFR current data

Sea surface currents consist of HFR-combined current data coming from two beam-forming WEllen RAdar stations operating in the Gulf of Trieste (Fig. 1b, red dots), the first located in Aurisina (Italy) and operated by the National Institute of Oceanography and Applied Geophysics (OGS) and the other in Piran (Slovenia), first operated by the National Institute of Biology (NIB) and later by the Slovenian Environmental Agency (ARSO). HFR works at a frequency of 24.5 MHz, with a space resolution of 1.5 km and a time resolution of 30 min. The data are freely available from the European HFR node at the following link: https://thredds.hfrnode.eu:8443/thredds/NRTcurrent/HFR-NAdr/HFR-NAdr_catalog.html (last access: 1 February 2023).

The quality-control standards from the EU high-frequency node (Corgnati et al., 2018) were applied to the data set. In addition, any remaining spikes were removed.

The grid points close to the HFR baseline have a major data coverage but show a more anisotropic behaviour than those far away (Fig. 1b). The P grid point was selected as the best candidate for the analysis of this paper, because it has a good data coverage, and the current increments are nearly isotropic.

2.2 Atmosphere data

The atmosphere data consist of the forecasted wind velocity field at 10 m above the surface in the Gulf of Trieste. These data are the output of the Weather Research and Forecasting (WRF) model (https://www.mmm.ucar.edu/models/wrf, last access: 29 March 2023), version 4.2.1, performed daily by the Agenzia Regionale per la Protezione dell'Ambiente del Friuli Venezia Giulia (ARPA FVG) using initial and boundary conditions from the National Oceanic and Atmospheric Administration Global Forecasting System (https://www.ncei.noaa.gov/products/weather-climate-models/global-forecast, last access: 29 March 2023). There are three calculation domains: the smallest one has the highest resolution, achieved with the two-way nesting technique, and provides the data used here (Goglio, 2018). The field has a spatial resolution of 2 km and a temporal resolution of 1 h.

Only the data from the closest WRF grid point to the P point are used for the following analysis. This point will be called “WRF_P” and it is shown in Fig. 1b.

3 Method: superstatistics

The analysis follows closely the framework introduced by Beck et al. (2005b). The reader familiar with it can skip the present section.

Consider a non-equilibrium system, characterized by a variable v and driven by an intensive parameter σ² (the temperature in Boltzmann statistics). This parameter σ² is approximately constant on a timescale T over which the system reaches a local equilibrium. The system fluctuates on a quicker decorrelation timescale τ≪T. Over long times, the shape of the probability distribution of v is non-trivial and depends on the probability distribution of σ² and on the type of interaction between σ² and v. The simplest example is the case of a Brownian particle moving in one dimension with a velocity v, a variance σ² which is related to the temperature and where the local conditional distribution p(v|σ²) is Gaussian (Maxwell–Boltzmann statistics in the canonical ensemble of statistical mechanics). According to the original point of view of Jaynes (1957), we interpret the local Gaussianity of the velocity, given a certain temperature or variance, as the least biased PDF estimate, as it is based on entropy maximization. Interpreting Shannon entropy as a measure of the lack of knowledge, it is possible to construct PDFs that maximize this lack of knowledge when some information is available. These obtained PDFs are the best estimates based on available information. We refer the reader to Sect. S5 in the Supplement to see how a continuous real variable with a given mean and variance maximizes the entropy if its PDF is Gaussian.

Figure 2 Classification of the wind regimes. (a) Daily wind (dots; the radial axes represent the speed and the azimuthal axes represent the direction) with the threshold wind speed (3 m s⁻¹, green line) and angle range (mistral in light blue, bora in blue, and sirocco in red) definitions. The bora has the highest speed peaks. (b) Wind regime type normalized histogram. The main strong wind regimes are bora and sirocco.

Superstatistics considers dynamics where the temperature or the variance is slowly varying in time. When the decorrelation time is much faster than the time over which the variance evolves (τ≪T), the dynamics results in a superposition of statistically different dynamics, so that the PDF of the velocity is

$$\begin{array}{}\text{(1)}& p\left(v\right)=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}f\left({\mathit{\sigma }}^{\mathrm{2}}\right)p\left(v\right|{\mathit{\sigma }}^{\mathrm{2}}\left)\mathrm{d}\right({\mathit{\sigma }}^{\mathrm{2}}),\end{array}$$

where

$$\begin{array}{}\text{(2)}& p\left(v\right|{\mathit{\sigma }}^{\mathrm{2}})={\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }{\mathit{\sigma }}^{\mathrm{2}}}}\right)}^{\mathrm{1}/\mathrm{2}}{e}^{-\frac{v^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}}\end{array}$$

is the local Gaussian PDF of the velocity v for a given σ² (the driving intensive parameter over T) and where f(σ²) is the PDF of σ².

Several studies (Beck, 2004; Beck et al., 2005a; Beck et al., 2005b) have shown that the thermodynamic beta $\mathit{\beta }=\mathrm{1}/{\mathit{\sigma }}^{\mathrm{2}}$ PDF often falls into three main classes: (i) χ²; (ii) inverse-χ²; (iii) log-normal. The justification for the classes is often scant, and other PDF shapes are also observed (Rizzo et al., 2004; Yalcin and Beck, 2013; Schäfer et al., 2021).

Starting from the longitudinal component u(t) and the latitudinal component v(t) of the sea surface currents in the P grid point in the Gulf of Trieste (Fig. 1b), we calculated the velocity increments $\mathit{\delta }u\left(t\right)=u(t+\mathit{\delta })-u\left(t\right)$ and $\mathit{\delta }v\left(t\right)=v(t+\mathit{\delta })-v\left(t\right)$, with the time increment $\mathit{\delta }={\mathrm{2}}^j\times \mathrm{30}\min$, $j=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{\dots },\mathrm{7}$, since the focus is on fluctuations and not on averages. The first objective is to extract, from each of these time series, the two superstatistical timescales τ_⋆ and T_⋆ (where ⋆ stands for δu or δv) with ${\mathit{\tau }}_{\star }<T_{\star }$ and to obtain ${\mathit{\sigma }}_{\star }^{\mathrm{2}}\left(t\right)$ and look at its PDF $f\left({\mathit{\sigma }}_{\star }^{\mathrm{2}}\right)$.

First, consider the timescale τ: the relaxation time τ_⋆ is defined as the exponential decay time of the autocorrelation function of the time series δu(t) (or δv(t)). For each δ, it is possible to calculate the autocorrelation function $C_{\mathit{\delta }u}\left(\mathrm{\Delta }t\right)=\langle \mathit{\delta }u(t^{\prime }+\mathrm{\Delta }t)\mathit{\delta }u\left(t^{\prime }\right)\rangle =\frac{{\int }_{\mathrm{0}}^{t_{max}-\mathrm{\Delta }t}\mathit{\delta }u(t^{\prime }+\mathrm{\Delta }t)\mathit{\delta }u\left(t^{\prime }\right)\mathrm{d}{t}^{\prime }}{{\int }_{\mathrm{0}}^{t_{max}-\mathrm{\Delta }t}\mathit{\delta }u\left(t^{\prime }\right)\mathit{\delta }u\left(t^{\prime }\right)\mathrm{d}{t}^{\prime }}$ (analogously for δv) and, from it, to extract the time for which $C_{\star }\left({\mathit{\tau }}_{\star }\right)=\mathrm{1}/e$. The time τ_⋆ is the decorrelation time of the signal.

Secondly, consider the timescale T: this timescale is the time for which the system is locally Gaussian. It is possible to calculate it through the average kurtosis value of the velocity increment as a function of time

$$\begin{array}{}\text{(3)}& {\mathit{\kappa }}_{\mathit{\delta }u}\left(\mathrm{\Delta }t\right)={\displaystyle \frac{\mathrm{1}}{t_{max}-\mathrm{\Delta }t}}\underset{\mathrm{0}}{\overset{t_{max}-\mathrm{\Delta }t}{\int }}{\displaystyle \frac{\langle {\left(\mathit{\delta }u-\stackrel{\mathrm{‾}}{\mathit{\delta }u}\right)}^{\mathrm{4}}{\rangle }_{t,\mathrm{\Delta }t}}{\langle {\left(\mathit{\delta }u-\stackrel{\mathrm{‾}}{\mathit{\delta }u}\right)}^{\mathrm{2}}{\rangle }_{t,\mathrm{\Delta }t}^{\mathrm{2}}}}\mathrm{d}t.\end{array}$$

Analogously for δv, where $\langle \mathrm{\cdots }{\rangle }_{t,\mathrm{\Delta }t}=\frac{\mathrm{1}}{\mathrm{\Delta }t}{\int }_t^{t+\mathrm{\Delta }t}\mathrm{\cdots }\mathrm{d}t$ is the internal Δt time average and $\stackrel{\mathrm{‾}}{\mathrm{\cdots }}$ is the total time average. By definition, T_⋆ is the time for which ${\mathit{\kappa }}_{\star }\left(T_{\star }\right)=\mathrm{3}$, which is the Gaussian value.

A large time separation ${\mathit{\tau }}_{\star }\ll T_{\star }$ is fundamental, because it guarantees that a local Gaussian equilibrium is reached. To determine the variance PDF $f\left({\mathit{\sigma }}_{\star }^{\mathrm{2}}\right)$, we also need a record length $T_{\text{obs}}\gg T_{\star }$.

As the reader can see from Eq. (2), the parameter ${\mathit{\sigma }}_{\star }^{\mathrm{2}}$ is precisely the variance of the local Gaussian distribution and, knowing T_⋆, it is possible to calculate ${\mathit{\sigma }}_{\star }^{\mathrm{2}}\left(t\right)$ as follows,

$$\begin{array}{}\text{(4)}& {\mathit{\sigma }}_{\mathit{\delta }u}^{\mathrm{2}}\left(t\right)=\langle \mathit{\delta }{u}^{\mathrm{2}}{\rangle }_{t,T}-\langle \mathit{\delta }u{\rangle }_{t,T}^{\mathrm{2}},\end{array}$$

and analogously for δv.

Superstatistics allows us to uncover the physics of a non-equilibrium system by finding the PDF of a variable by separating the timescales and bringing out the evolution of the local equilibrium PDF.

4 Results and discussion

First we have identified and clustered the wind regimes blowing over the Gulf of Trieste in the analysed period. Since bora and sirocco are synoptic winds, the daily wind time series have been used (Fig. 2a). For each day, a wind regime is attributed: (i) low wind if the wind speed is < 3 m s⁻¹ (green area); (ii) mistral if the wind speed is ≥ 3 m s⁻¹ and its direction θ is −67.5^∘ N ≤θ  < 22.5^∘ N (light-blue area); (iii) bora if the wind speed is ≥ 3 m s⁻¹ and its direction θ is 22.5^∘ N ≤θ < 112.5^∘ N (blue area); and (iv) sirocco if the wind speed is ≥ 3 m s⁻¹ and its direction θ is 112.5^∘ N $\le \mathit{\theta }\le$ 206^∘ N (red area).

It is interesting to see that the bora shows the highest speed peaks (Fig. 2a), re-veiling the strongest wind forcing over the Gulf of Trieste, with a maximum daily wind speed of about 50 km h⁻¹. About the wind regime occurrence, as can be seen from the histogram in Fig. 2b, the low wind regime is the most frequent (occurrence of more than 60 %), while among strong wind regimes the bora arises more often (occurrence of more than 25 %). The sirocco develops less frequently (occurrence of just below 10 %), while the mistral has an occurrence of less than 1 % (it counts just five daily events in a period of almost 2 years; Fig. 2a), providing insufficient statistics, so it will be ignored in the rest of the analysis.

Next, we have considered the current data. Before following the τ_⋆ and T_⋆ calculation procedures, the velocity increments PDFs were calculated: all these PDFs are non-Gaussian and fat-tailed with kurtosis values greater than 3. The kurtosis value decreases with an increasing δ, approaching the value 3, showing the fact that the larger the increment time, the more the uncorrelated velocity increments are and the closer the time series is to a Gaussian process.

Figure 3 Velocity increment kurtosis value κ_⋆ as a function of time Δt from the P grid point and T_⋆ estimation for the δu (a) and δv (b) velocity increments (δ and T_⋆ in units of 30 min). The timescale T_⋆ always takes values longer than 10 h, showing the timescale separation as τ_⋆ assumes values between 15 min and 3 h.

The velocity increments correlation functions were calculated together with the relaxation time τ_⋆ estimations, as explained in Sect. 3. The relaxation time τ_⋆ assumes values between approximately 15 min and 3 h, increasing with an increasing δ. The velocity increments kurtosis values in function of the time were calculated together with the Gaussianity timescale T_⋆ estimations (Fig. 3), as explained in Sect. 3. The timescale T_⋆ takes values longer than 10 h, increasing with δ. This fact shows that the larger the increment time, the more uncorrelated the time series is. We study the horizontal turbulent dynamics of the upper sea layer, evolving on timescales longer than a few hours, subject to the synoptic atmospheric forcing, with a typical timescale of several days. Choosing the velocity time increment $\mathit{\delta }=\mathrm{8}\times$ 30 min = 4 h leads to $T_{\star }\simeq$ 2 d, which is in agreement with the ocean surface-layer turbulent and the synoptic atmospheric timescales, respectively. In the following all the results will refer to this velocity time increment.

Figure 4 Velocity increments δu(t), δv(t) (in blue), and their respective ${\mathit{\sigma }}_{\star }^{\mathrm{2}}\left(t\right)$ (in orange) time series with $\mathit{\delta }=\mathrm{8}\times$ 30 min = 4 h from the P grid point and wind regime periods from the WRF_P grid point (coloured shades). The vertical scales are different in the two plots for better visualization.

The resultant ${\mathit{\sigma }}_{\star }^{\mathrm{2}}\left(t\right)$ time series was calculated, according to Eq. (4). It is possible to see that ${\mathit{\sigma }}_{\star }^{\mathrm{2}}\left(t\right)$ varies more slowly than the velocity increment (Fig. 4). This fact is also apparent in the comparison between the velocity increments autocorrelation and the ${\mathit{\sigma }}_{\star }^{\mathrm{2}}$ autocorrelation (not shown). It is possible to see that ${\mathit{\sigma }}_{\mathit{\delta }u}^{\mathrm{2}}$ has higher peaks with respect to ${\mathit{\sigma }}_{\mathit{\delta }v}^{\mathrm{2}}$. This means that δu locally has a higher variance and thus a greater variability than δv. This is probably due to the geometry of the Gulf of Trieste whose western side is open. However this difference is not particularly pronounced, in fact it does not exceed 1 order of magnitude. Moreover no clear seasonal cycle is evident.

From our estimations of the observations (not shown), the thermodynamic betas ${\mathit{\beta }}_{\star }=\mathrm{1}/{\mathit{\sigma }}_{\star }^{\mathrm{2}}$ do not fall into the χ² and log-normal universality classes, but the fits suggest an inverse-Gamma distribution (a universality class proposed by Beck et al., 2005b). In general, the inverse-Gamma distribution with the shape parameter n has divergent nth and higher-order moments, bringing limitations in fitting data for low n. As we will discuss below, we identify the shape parameter n with the degrees of freedom of the system, in our case equal to 2, leading to an analytical inverse-Gamma distribution for β_⋆ with only the first moment defined. For this reason we propose considering ${\mathit{\sigma }}_{\star }^{\mathrm{2}}={\mathit{\beta }}_{\star }^{-\mathrm{1}}$ instead of β_⋆. The PDF of the variance $f\left({\mathit{\sigma }}_{\star }^{\mathrm{2}}\right)$ is a Gamma distribution (Sect. S1) with a fixed shape parameter equal to 2, where all positive moments converge:

$$\begin{array}{}\text{(5)}& f\left({\mathit{\sigma }}_{\star }^{\mathrm{2}}\right)={\mathrm{\Gamma }}_{\mathrm{2},{\mathit{\lambda }}_{\star }}\left({\mathit{\sigma }}_{\star }^{\mathrm{2}}\right)={\mathit{\lambda }}_{\star }^{\mathrm{2}}{\mathit{\sigma }}_{\star }^{\mathrm{2}}{e}^{-{\mathit{\lambda }}_{\star }{\mathit{\sigma }}_{\star }^{\mathrm{2}}}.\end{array}$$

Our environmental system has 2 spatial degrees of freedom, so we can describe ${\mathit{\sigma }}_{\star }^{\mathrm{2}}$ as the sum of two independent positive-defined variables (A_⋆ and B_⋆). Again, we apply entropy maximization. Each of these independent variables is distributed to maximize entropy: since it is a positive variable, it is distributed as an exponential distribution, a particular case of the Gamma distribution (Sect. S5):

$$\begin{array}{}\text{(6)}& \begin{array}{rl}& {\mathit{\sigma }}_{\star }^{\mathrm{2}}=\underset{A_{\star }\text{ and }{B}_{\star }\text{ are independent}}{\underbrace{A_{\star }+B_{\star };}}\\ & \text{ where }\left\{\begin{array}{l}{A}_{\star }\sim {\mathrm{\Gamma }}_{\mathrm{1},{\mathit{\lambda }}_{\star }}\left(x\right)={\mathit{\lambda }}_{\star }{e}^{-{\mathit{\lambda }}_{\star }x}\\ B_{\star }\sim {\mathrm{\Gamma }}_{\mathrm{1},k_{\star }}\left(x\right)=k_{\star }{e}^{-k_{\star }x}.\end{array}\right.\end{array}\end{array}$$

This fact would give the following ${\mathit{\sigma }}_{\star }^{\mathrm{2}}$ distribution (Sect. S2):

$$\begin{array}{}\text{(7)}& f\left({\mathit{\sigma }}_{\star }^{\mathrm{2}}\right)={\displaystyle \frac{{\mathit{\lambda }}_{\star }{k}_{\star }}{{\mathit{\lambda }}_{\star }-k_{\star }}}(e^{-k_{\star }{\mathit{\sigma }}_{\star }^{\mathrm{2}}}-e^{-{\mathit{\lambda }}_{\star }{\mathit{\sigma }}_{\star }^{\mathrm{2}}}).\end{array}$$

However, since the data are nearly isotropic (Fig. 1b), we can suppose the equality of the coefficient for both degrees of freedom ${\mathit{\lambda }}_{\star }=k_{\star }$ to obtain Eq. (5).

Figure 5 Velocity increment PDFs (δu on (a, c) and δv on (b, d) from the P grid point and $\mathit{\delta }=\mathrm{8}\times$ 30 min with the best fit from Eq. (8). The PDFs are reported in a lin–lin plot (a, b) to enhance the peaks and in a log–log plot(c, d) to enhance the tails. The thin spotted line represents the HFR data, the thick line represents the fit, and the different colours are used to discriminate between the different wind regimes.

Using Eqs. (1), (2), and (5), it is possible to calculate the PDF of the original signal (Sects. S3 and S4)

$$\begin{array}{}\text{(8)}& \begin{array}{rl}p\left(\mathit{\delta }u\right)=& \underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}f\left({\mathit{\sigma }}_{\mathit{\delta }u}^{\mathrm{2}}\right)p\left(\mathit{\delta }u|{\mathit{\sigma }}_{\mathit{\delta }u}^{\mathrm{2}}\right)\mathrm{d}\left({\mathit{\sigma }}_{\mathit{\delta }u}^{\mathrm{2}}\right)\\ =& {\displaystyle \frac{\sqrt{\mathrm{2}{\mathit{\lambda }}_{\mathit{\delta }u}}{e}^{-\sqrt{\mathrm{2}{\mathit{\lambda }}_{\mathit{\delta }u}}\left|\mathit{\delta }u\right|}\left(\sqrt{\mathrm{2}{\mathit{\lambda }}_{\mathit{\delta }u}}\right|\mathit{\delta }u|+\mathrm{1})}{\mathrm{4}}}\end{array}\end{array}$$

and analogously for δv. We emphasize that this PDF is obtained by maximizing entropy twice: for the Gaussian-distributed $p\left(\mathit{\delta }u\right|{\mathit{\sigma }}_{\mathit{\delta }u}^{\mathrm{2}})$ and for each degree of freedom of $f\left({\mathit{\sigma }}_{\star }^{\mathrm{2}}\right)$. We have fitted Eq. (8) on the observed p(δu) and p(δv) and their PDFs during the different wind regimes separately, obtaining for each curve the best coefficient λ_⋆ (Fig. 5). As can be seen in Fig. 5, the analytical PDF expression is valid for the entire velocity increment data set as well as for the different wind regimes separately. From the expression in Eq. (8), it is possible to also derive the analytical expression of the velocity increment second-order moments $s_{\star }^{\mathrm{2}}$ (Sects. S3 and S4):

$$\begin{array}{}\text{(9)}& \left.\begin{array}{c}\text{var}\left(\mathit{\delta }u\right)\\ \text{var}\left(\mathit{\delta }v\right)\end{array}\right\}=s_{\star }^{\mathrm{2}}={\displaystyle \frac{\mathrm{2}}{{\mathit{\lambda }}_{\star }}}=\sqrt{\mathrm{2}}{\mathit{\sigma }}_{\star \mathrm{\Gamma }}=\mathrm{2}{\mathit{\sigma }}_{\star \text{exp}}.\end{array}$$

which are analytically related to the ${\mathrm{\Gamma }}_{\mathrm{2},{\mathit{\lambda }}_{\star }}\left({\mathit{\sigma }}_{\star }^{\mathrm{2}}\right)$ standard deviation (σ_⋆Γ) and to the exponential distribution standard deviation (σ_⋆exp). The velocity increment variance $s_{\star }^{\mathrm{2}}$ is twice the exponential distribution standard deviation σ_⋆exp, as there are 2 spatial degrees of freedom. Equation (8) can then be expressed as

$$\begin{array}{}\text{(10)}& p\left(\mathit{\delta }u\right)={\displaystyle \frac{e^{-\mathrm{2}\left|\mathit{\delta }u\right|/s_{\mathit{\delta }u}}\left(\frac{\mathrm{2}\left|\mathit{\delta }u\right|}{s_{\mathit{\delta }u}}+\mathrm{1}\right)}{\mathrm{2}{s}_{\mathit{\delta }u}}}.\end{array}$$

The s_⋆ standard deviation numerical values from the best-fit coefficients λ_⋆ are reported in Table 1.

Table 1 Best-fit standard deviations from Eq. (9). The bora, the strongest wind forcing, leads to the lowest fluctuations of the velocity increments.

Figure 6 PDFs of ${\mathit{\sigma }}_{\mathit{\delta }u}^{\mathrm{2}}$ (a, c) and ${\mathit{\sigma }}_{\mathit{\delta }v}^{\mathrm{2}}$ (b, d) with $\mathit{\delta }=\mathrm{8}\times$ 30 min at the P HFR grid point: total $f\left({\mathit{\sigma }}_{\star }^{\mathrm{2}}\right)$ and clustered $f\left({\mathit{\sigma }}_{\star }^{\mathrm{2}}\right|e)$ by the wind regime e. For the clustering operation, the ${\mathit{\sigma }}_{\star }^{\mathrm{2}}$ time series were shifted by $T_{\star }/\mathrm{2}$, and then the values were collected as a function of the daily wind regime. The PDFs are reported in a lin–lin plot (a, b) and in a log–log plot (c, d). The thin spotted line represents the HFR data, the thick line represents the analytical distribution, and the different colours are used to discriminate between the different wind regimes.

Using the best-fit coefficients λ_⋆ just derived, the solution of Eq. (5) is plotted and compared to the observed ${\mathit{\sigma }}_{\star }^{\mathrm{2}}$ PDFs in Fig. 6.

5 Conclusions

We have observed that the sea surface current velocity increments show fat tails in their distributions. Superstatistics is a powerful tool for determining the PDF. In practice, it relies on the timescale separation between the fast fluctuations (for which the Maxwell–Boltzmann statistics applies) and their slow driver (the temperature in statistical thermodynamics), which makes the local PDF evolve with time. From the velocity increment time series, we have extracted the two different timescales: the relaxation time ${\mathit{\tau }}_{\star }\simeq$ 1 h 50 min and the larger $T_{\star }\simeq$ 2 d, which is the Gaussianity timescale. According to the original point of view of Jaynes (1957), the local Gaussianity of the velocity increments is the least-biased PDF estimate based on the information available as it reflects the maximization of entropy. From the T_⋆ timescale we have extracted the local Gaussian variance time series ${\mathit{\sigma }}_{\star }^{\mathrm{2}}\left(t\right)$ (the “temperature”). The maximum-entropy PDF for a positive variable, as a variance, is exponential. In our case of 2 degrees of freedom we obtained the convolution of the two identical exponentials, i.e. a Gamma ${\mathrm{\Gamma }}_{\mathrm{2},{\mathit{\lambda }}_{\star }}$ distribution function. Note that Beck and Cohen (2003) interprets the degrees of freedom as the number of the elements in the sum of squared identically distributed Gaussian variables, which give a Gamma-distributed variable. The degrees of freedom in their interpretation correspond to twice our value. The ${\mathit{\sigma }}_{\star }^{\mathrm{2}}$ Gamma distribution corresponds to the inverse-Gamma β_⋆ class described by Beck et al. (2005b). Nevertheless, this distribution, representing a system with n degrees of freedom, presents divergent m≥n-order moments, bringing limitations in fitting data for low n, as in the present case. This problem completely disappears when ${\mathit{\sigma }}_{\star }^{\mathrm{2}}={\mathit{\beta }}_{\star }^{-\mathrm{1}}$ is considered.

From the PDF of the variance we obtained the analytical distribution of the velocity increments. Fitting the distribution to our data we observe a good agreement, even if superstatistics works in the limit of infinite timescale separations ${\mathit{\tau }}_{\star }\ll T_{\star }$ (a factor of few tens in our case) and $T_{\star }\ll T_{\text{obs}}$ (a factor of a few hundreds in our case). When the analysis is performed for the different wind regimes separately, the same analytical laws of the PDFs are obtained with different variance (Table 1), which is the only free parameter in the fit. The good agreement also for the different wind regimes is surprising as the length of the wind regimes is comparable to T_⋆ and for Sirocco winds we have only a few tens of independent events to fit the PDF of the variance. The good fits for the different wind regimes with the same analytic form points towards a universal behaviour. We showed that superstatistics is more than a tool for fitting fat-tailed PDFs: it gives the characteristic timescales and when combined with the maximum-entropy principle it indicates the number of the degrees of freedom.

The velocity increments second-order moments differ between the forcing regimes (bora, sirocco, and low wind), with the strongest wind forcing (bora) leading to the lowest fluctuations of the velocity increments. This might be a consequence of the sea surface current tendency to settle on a mean outflow (from the gulf) during the bora or might indicate that an increased wave activity and turbulent exchange in the vertical of the horizontal momentum, operating at even shorter timescales than δ, lead to a damping of the sea surface current fluctuations.

Parameterizations of turbulent motions generally rely on a velocity scale and a length scale, averaged over a large timescale. The availability of a considerable amount of data allows us here to go further by determining the PDF of the velocity fluctuations at different timescales. More precisely, the PDF of the fast velocity fluctuations is Gaussian and the PDF of the velocity variance at timescales larger than T is Γ_2,λ. As already mentioned, for the Gulf of Trieste T≃ 2 d for velocity increments over a few hours (Fig. 3). The analytic form of the combined PDF of the velocity increments (Eq. 8) can be used to develop parameterizations of the dynamics in the ocean surface layer at timescales faster than a few days and can discriminate between existing parameterizations. Our findings show that parameterizations of surface current fluctuations at timescales shorter than T∼ 2 d can rely on Gaussian fluctuations. For longer timescales the dynamics is non-Gaussian and fat-tailed, meaning that extreme events are more ponderous and a parameterization has to account for it. The findings currently guide us in constructing stochastic differential equations that will lead to the same PDFs obtained here from a superstatistical analysis of the data. When the stochastic differential equation is known, it is mathematically straightforward to derive parameterizations of the turbulent dynamics.