Articles | Volume 30, issue 4
Research article
21 Nov 2023
Research article |  | 21 Nov 2023

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

Sofia Flora, Laura Ursella, and Achim Wirth

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 β=1/σ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 δ=8×30min=4h 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 “WRFP”. 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−1 and bringing cold and dry continental air from the north-east over the sea (Poulain and Raicich2001). 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 Raicich2001). 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 (Crombie1955). 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: (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 (, 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 (, 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 (Goglio2018). 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 “WRFP” 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 σ2 (the temperature in Boltzmann statistics). This parameter σ2 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 σ2 and on the type of interaction between σ2 and v. The simplest example is the case of a Brownian particle moving in one dimension with a velocity v, a variance σ2 which is related to the temperature and where the local conditional distribution p(v|σ2) 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 2Classification 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−1, 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

(1) p ( v ) = 0 f ( σ 2 ) p ( v | σ 2 ) d ( σ 2 ) ,


(2) p ( v | σ 2 ) = 1 2 π σ 2 1 / 2 e - v 2 2 σ 2

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

Several studies (Beck, 2004; Beck et al., 2005a; Beck et al., 2005b) have shown that the thermodynamic beta β=1/σ2 PDF often falls into three main classes: (i) χ2; (ii) inverse-χ2; (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 δu(t)=u(t+δ)-u(t) and δv(t)=v(t+δ)-v(t), with the time increment δ=2j×30min, j=0,1,2,,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 τ<T and to obtain σ2(t) and look at its PDF f(σ2).

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δu(Δt)=δu(t+Δt)δu(t)=0tmax-Δtδu(t+Δt)δu(t)dt0tmax-Δtδu(t)δu(t)dt (analogously for δv) and, from it, to extract the time for which C(τ)=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

(3) κ δ u ( Δ t ) = 1 t max - Δ t 0 t max - Δ t δ u - δ u 4 t , Δ t δ u - δ u 2 t , Δ t 2 d t .

Analogously for δv, where t,Δt=1Δttt+Δtdt is the internal Δt time average and is the total time average. By definition, T is the time for which κ(T)=3, which is the Gaussian value.

A large time separation τT is fundamental, because it guarantees that a local Gaussian equilibrium is reached. To determine the variance PDF f(σ2), we also need a record length TobsT.

As the reader can see from Eq. (2), the parameter σ2 is precisely the variance of the local Gaussian distribution and, knowing T, it is possible to calculate σ2(t) as follows,

(4) σ δ u 2 ( t ) = δ u 2 t , T - δ u t , T 2 ,

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−1 (green area); (ii) mistral if the wind speed is  3 m s−1 and its direction θ is 67.5 N θ  < 22.5 N (light-blue area); (iii) bora if the wind speed is  3 m s−1 and its direction θ is 22.5 N θ< 112.5 N (blue area); and (iv) sirocco if the wind speed is  3 m s−1 and its direction θ is 112.5 N θ 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−1. 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 3Velocity 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 δ=8× 30 min= 4 h leads to T 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 4Velocity increments δu(t), δv(t) (in blue), and their respective σ2(t) (in orange) time series with δ=8× 30 min= 4 h from the P grid point and wind regime periods from the WRFP grid point (coloured shades). The vertical scales are different in the two plots for better visualization.


The resultant σ2(t) time series was calculated, according to Eq. (4). It is possible to see that σ2(t) varies more slowly than the velocity increment (Fig. 4). This fact is also apparent in the comparison between the velocity increments autocorrelation and the σ2 autocorrelation (not shown). It is possible to see that σδu2 has higher peaks with respect to σδv2. 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 β=1/σ2 do not fall into the χ2 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 σ2=β-1 instead of β. The PDF of the variance f(σ2) is a Gamma distribution (Sect. S1) with a fixed shape parameter equal to 2, where all positive moments converge:

(5) f σ 2 = Γ 2 , λ σ 2 = λ 2 σ 2 e - λ σ 2 .

Our environmental system has 2 spatial degrees of freedom, so we can describe σ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):

(6) σ 2 = A + B ; A  and  B  are independent  where  A Γ 1 , λ ( x ) = λ e - λ x B Γ 1 , k ( x ) = k e - k x .

This fact would give the following σ2 distribution (Sect. S2):

(7) f ( σ 2 ) = λ k λ - k ( e - k σ 2 - e - λ σ 2 ) .

However, since the data are nearly isotropic (Fig. 1b), we can suppose the equality of the coefficient for both degrees of freedom λ=k to obtain Eq. (5).

Figure 5Velocity increment PDFs (δu on (a, c) and δv on (b, d) from the P grid point and δ=8× 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)

(8) p ( δ u ) = 0 f σ δ u 2 p δ u | σ δ u 2 d σ δ u 2 = 2 λ δ u e - 2 λ δ u | δ u | ( 2 λ δ u | δ u | + 1 ) 4

and analogously for δv. We emphasize that this PDF is obtained by maximizing entropy twice: for the Gaussian-distributed p(δu|σδu2) and for each degree of freedom of f(σ2). 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 s2 (Sects. S3 and S4):

(9) var ( δ u ) var ( δ v ) = s 2 = 2 λ = 2 σ Γ = 2 σ exp .

which are analytically related to the Γ2,λ(σ2) standard deviation (σ⋆Γ) and to the exponential distribution standard deviation (σ⋆exp). The velocity increment variance s2 is twice the exponential distribution standard deviation σ⋆exp, as there are 2 spatial degrees of freedom. Equation (8) can then be expressed as

(10) p ( δ u ) = e - 2 | δ u | / s δ u 2 | δ u | s δ u + 1 2 s δ u .

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

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

Download Print Version | Download XLSX

Figure 6PDFs of σδu2 (a, c) and σδv2 (b, d) with δ=8× 30 min at the P HFR grid point: total f(σ2) and clustered f(σ2|e) by the wind regime e. For the clustering operation, the σ2 time series were shifted by T/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 σ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 τ 1 h 50 min and the larger T 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 σ2(t) (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 Γ2,λ 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 σ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 mn-order moments, bringing limitations in fitting data for low n, as in the present case. This problem completely disappears when σ2=β-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 τT (a factor of few tens in our case) and TTobs (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.

Data availability

The HFR sea surface current data of the Gulf of Trieste are freely available from the European HFR node at the following DOI: (Ursella et al.2023). The WRF forecasted wind field is obtainable upon request from ARPA FVG (, CRMA2023).


The supplement related to this article is available online at:

Author contributions

SF performed the major part of the research and the writing of the manuscript. LU and AW contributed to the research and the writing.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.


The forecasts, analyses, and related services are based on data and products of the Regional Center for Environmental Modelling (CRMA), which is a sector of the Friuli Venezia Giulia Environmental Agency (ARPA FVG) ITALY (, last access: 29 March 2023). The Piran HFR station has been operated by the Slovenian Environment Agency since July 2022. Between 2015 and July 2022 it was operated by the National Institute of Biology. Sofia Flora acknowledges the hospitality of the LEGI laboratory (Grenoble, France). The authors would like to thank the anonymous reviewers who significantly improved the quality of the paper.

Financial support

This work was funded by the Initiatives de recherche À Grenoble Alpes (IRGA) 2022 through the FASIOM project.

Review statement

This paper was edited by Norbert Marwan and reviewed by three anonymous referees.


Beck, C.: Superstatistics: theory and applications, Continuum Mech. Therm., 16, 293–304, 2004. a, b, c

Beck, C. and Cohen, E. G.: Superstatistics, Physica A, 322, 267–275, 2003. a, b

Beck, C., Cohen, E., and Rizzo, S.: Atmospheric turbulence and superstatistics, Europhysics News, 36, 189–191, 2005a. a, b

Beck, C., Cohen, E. G., and Swinney, H. L.: From time series to superstatistics, Phys. Rev. E, 72, 056133, 2005b. a, b, c, d, e, f, g, h

Bogunović, B. and Malačič, V.: Circulation in the Gulf of Trieste: Measurements and model results, Nuovo Cimento C, 31, 301–326, 2008. a

CRMA – Centro Regionale di Modellistica Ambientale: Progetto NAUSICA – Downscaling di analisi meteorologiche ad alta risoluzione sul dominio Alpe Adria,, last access: 29 March 2023. a

Corgnati, L., Mantovani, C., Novellino, A., Rubio, A., and Mader, J.: Recommendation Report 2 on improved common procedures for HFR QC analysis, JERICO-NEXT WP5-Data Management, Deliverable 5.14, Version 1.0., FREMER,, 2018. a

Cosoli, S., Gačić, M., and Mazzoldi, A.: Surface current variability and wind influence in the northeastern Adriatic Sea as observed from high-frequency (HF) radar measurements, Cont. Shelf Res., 33, 1–13, 2012. a, b, c, d, e

Cosoli, S., Ličer, M., Vodopivec, M., and Malačič, V.: Surface circulation in the Gulf of Trieste (northern Adriatic Sea) from radar, model, and ADCP comparisons, J. Geophys. Res.-Oceans, 118, 6183–6200, 2013. a, b

Costa, M. O., Silva, R., and Anselmo, D. H. A. L.: Superstatistical and DNA sequence coding of the human genome, Phys. Rev. E, 106, 064407,, 2022. a

Crombie, D. D.: Doppler spectrum of sea echo at 13.56 Mc./s., Nature, 175, 681–682, 1955. a

dos Santos, M., Menon, L., and Cius, D.: Superstatistical approach of the anomalous exponent for scaled Brownian motion, Chaos Soliton. Fract., 164, 112740,, 2022. a

Ghil, M., Yiou, P., Hallegatte, S., Malamud, B. D., Naveau, P., Soloviev, A., Friederichs, P., Keilis-Borok, V., Kondrashov, D., Kossobokov, V., Mestre, O., Nicolis, C., Rust, H. W., Shebalin, P., Vrac, M., Witt, A., and Zaliapin, I.: Extreme events: dynamics, statistics and prediction, Nonlin. Processes Geophys., 18, 295–350,, 2011. a, b

Goglio, A. C.: Progetto NAUSICA – Downscaling di analisi meteorologiche ad alta risoluzione sul dominio Alpe Adria, ARPA FVG, (last access: 29 March 2023), 2018. a

Jaynes, E. T.: Information theory and statistical mechanics, Phys. Rev., 106, 620–630, 1957. a, b, c

Jeffries, M. A. and Lee, C. M.: A climatology of the northern Adriatic Sea's response to bora and river forcing, J. Geophys. Res.-Oceans, 112, C03S02,, 2007. a, b

Lorente, P., Aguiar, E., Bendoni, M., Berta, M., Brandini, C., Cáceres-Euse, A., Capodici, F., Cianelli, D., Ciraolo, G., Corgnati, L., Dadić, V., Doronzo, B., Drago, A., Dumas, D., Falco, P., Fattorini, M., Gauci, A., Gómez, R., Griffa, A., Guérin, C.-A., Hernández-Carrasco, I., Hernández-Lasheras, J., Ličer, M., Magaldi, M. G., Mantovani, C., Mihanović, H., Molcard, A., Mourre, B., Orfila, A., Révelard, A., Reyes, E., Sánchez, J., Saviano, S., Sciascia, R., Taddei, S., Tintoré, J., Toledo, Y., Ursella, L., Uttieri, M., Vilibić, I., Zambianchi, E., and Cardin, V.: Coastal high-frequency radars in the Mediterranean – Part 1: Status of operations and a framework for future development, Ocean Sci., 18, 761–795,, 2022. a

Malačič, V., Petelin, B., Gačić, M., Artegiani, A., and Orlić, M.: Regional Studies, Springer Netherlands, Dordrecht,, pp. 167–216, 2001. a, b

Okorie, U. S., Ikot, A. N., Okon, I. B., Obagboye, L. F., Horchani, R., Abdullah, H. Y., Qadir, K. W., and Abdel-Aty, A.-H.: Exact solutions of κ-dependent Schrödinger equation with quantum pseudo-harmonic oscillator and its applications for the thermodynamic properties in normal and superstatistics, Sci. Rep.-UK, 13, 2108,, 2023. a

Poulain, P.-M. and Raicich, F.: Forcings, Springer Netherlands, Dordrecht,, pp. 45–65, 2001. a, b, c, d

Pullen, J., Doyle, J. D., Haack, T., Dorman, C., Signell, R. P., and Lee, C. M.: Bora event variability and the role of air-sea feedback, J. Geophys. Res.-Oceans, 112, C03S18, 2007. a, b

Querin, S., Crise, A., Deponte, D., and Solidoro, C.: Numerical study of the role of wind forcing and freshwater buoyancy input on the circulation in a shallow embayment (Gulf of Trieste, Northern Adriatic Sea), J. Geophys. Res.-Oceans, 111, C03S16,, 2006. a, b, c, d

Querin, S., Cosoli, S., Gerin, R., Laurent, C., Malačič, V., Pristov, N., and Poulain, P.-M.: Multi-platform, high-resolution study of a complex coastal system: The TOSCA experiment in the Gulf of Trieste, J. Mar. Sci. Eng. 9, 5,, 2021. a, b, c, d

Raicich, F., Malačič, V., Celio, M., Giaiotti, D., Cantoni, C., Colucci, R., Čermelj, B., and Pucillo, A.: Extreme air-sea interactions in the Gulf of Trieste (North Adriatic) during the strong Bora event in winter 2012, J. Geophys. Res.-Oceans, 118, 5238–5250, 2013. a, b

Rizzo, S., Rapisarda, A., and Group, C.: Environmental atmospheric turbulence at Florence airport, in: AIP Conference Proceedings, vol. 742, 176–181, American Institute of Physics,, 2004. a, b

Schäfer, B., Heppell, C. M., Rhys, H., and Beck, C.: Fluctuations of water quality time series in rivers follow superstatistics, Iscience, 24, 102881,, 2021. a, b, c

Ursella, L., Čermelj, B., Ličer, M., and Martini, S.: HFR-NAdr (High Frequency Radar NAdr network), European HFRadar Node [data set],, 2023. a

Weber, J., Reyers, M., Beck, C., Timme, M., Pinto, J. G., Witthaut, D., and Schäfer, B.: Wind power persistence characterized by superstatistics, Sci. Rep.-UK, 9, 1–15, 2019.  a

Yalcin, G. C. and Beck, C.: Environmental superstatistics, Physica A, 392, 5431–5452, 2013. a, b

Short summary
An increasing amount of data allows us to move from low-order moments of fluctuating observations to their PDFs. We found the analytical fat-tailed PDF form (a combination of Gaussian and two-exponential convolutions) for 2 years of sea surface current increments in the Gulf of Trieste, using superstatistics and the maximum-entropy principle twice: on short and longer timescales. The data from different wind regimes follow the same analytical PDF, pointing towards a universal behaviour.