Estimating the population mean for a vertical profile of energy dissipation rate.

Energy dissipation rates are an important characteristic of turbulence; however, their magnitude in observational profiles can be incorrectly determined owing to their irregular appearance during vertical evolution. By analysing the data obtained from oceanic turbulence measurements, we demonstrate that the vertical sequences of energy dissipation rates exhibit a scaling property. Utilising this property, we propose a method to estimate the population mean for a profile. For scaling in the observed profiles, we demonstrate that our data exhibit a statistical property consistent with that exhibited by the universal multifractal model. Meanwhile, the population mean and its uncertainty can be estimated by inverting the probability distribution obtained by Monte Carlo simulations of a cascade model; to this end, observational constraints from several moments are imposed over each vertical sequence. This approach enables us to determine, to some extent, whether a profile shows an occasionally large mean or whether the population mean itself is large. Thus, it will contribute to the refinement of the regional estimation of the ocean energy budget, where only a small amount of turbulence observation data is available.

Introduction

Numerous existing studies have highlighted the importance of determining energy dissipation rates to investigate ocean general circulation[1,2]. Therefore, several observational studies have been conducted to obtain the vertical profiles of energy dissipation rates using ocean microstructure profilers[3,4]. In addition, to understand the statistics of the irregular evolution of observational profiles, studies have been conducted from the perspective of statistical fluid mechanics, as summarised below.

In fully developed turbulence, an inertial subrange of length scales exists wherein the advective term dominates the molecular viscosity term in the Navier–Stokes equation[5]. In this inertial subrange, a cascade of energies can be observed from large to small scales, as intuitively stated by Richardson[6]. In the first quantitative theory on energy cascades, Kolmogorov[7] established a relationship wherein velocity fluctuations are locally isotropic and are determined by the homogeneous energy dissipation rate; here, homogeneous means that the statistical property is independent of the position x,where v denotes the velocity; , the energy dissipation rate; , the distance between the points; and , the expected value. Subsequently, the energy dissipation rate was argued to vary, exhibiting considerable random fluctuations[8]. Thus, a refined theory[9] was proposed to address this issue. This theory stated that (1) , which is the logarithm of the spatially averaged energy dissipation rate over scale r, obeys a Gaussian distribution, and (2) its variance obeys , where L denotes the outer scale; A, a constant associated with the macrostructure of flow; and , the intermittency constant.

In addition, several experimental studies[10,11] demonstrated that small-scale dissipation is a random field that has a spatial structure with power–law correlations,Then, Yaglom[12] formulated a quantitative model, which was consistent with the log-normal scaling presented by Kolmogorov[9] and the power–law correlations, as a multiplicative cascade, where was expressed with a binary tree comprising independent and identically distributed (i.i.d) random variables, (),where are the positions with equal spacing and is the floor function, which assigns the integer that satisfies . If the random variables are set to have the moment exponent , then the energy conservation in a probabilistic sense, , and the log-normal scaling in Kolmogorov[9] are reproduced. Moreover, correlation (2) is reproduced because we have , where for small r[12,13].

Several alternative multiplicative cascade models have been developed with different generators, including the model[14], random model[15], model[16], p model[17], log-stable model[18], and log-Poisson model[19]. An important observation regarding Yaglom’s cascade is that the property required for the law of random variable W can be formulated such that the product of several random variables still obeys the same class of distribution, , with [20]. Consistent with this condition, the universal multifractal model[18] employs a stable Lévy generator, , that is maximally left skewed and satisfies . This results in a simple and nonanalytic form of the moment exponent, where is the multifractal index, which can be a non-integer, and is the codimension of the mean. The universal multifractal model is the most promising model. This model can well reproduce the variability in several phenomena including turbulence, other geophysical phenomena, and several fractal-like appearances in natural and man-made objects.

Based on this theory, we discuss a refined statistical treatment of the vertical profiles of the observed energy dissipation rates. We first distinguish the ’mean energy dissipation rate’, which refers to the sample (arithmetic) mean over a profile, and the ’energy input rate’, which refers to the population mean for a profile.

Thus, we reconsider one of the basic questions in the observational study of ocean turbulence: When a vertical profile of the energy dissipation rate is given, how can one estimate the energy input rate or the population mean of the energy dissipation rate for a profile, which has been commonly equated with the arithmetic mean over the profile? Our question pertains to whether one can obtain information regarding the energy input rate beyond the arithmetic mean. The answer is yes, because we can construct a model for the turbulent cascade process and solve the inversion problem to obtain the energy input rate under an observational constraint. In this study, we first show that the observed profiles of the depth-averaged energy dissipation rate, , exhibit a scaling property consistent with that of the universal multifractal model. Then, we construct a multiplicative cascade simulation model that describes the statistics of the observational data. Finally, we propose a method to explain certain statistics of the observed profiles based on the simulation model and develop an inversion method to estimate the energy input rate. This result illustrates a systematic method of gaining further quantitative information from profile data.

Methods

Observational data

In this section, we describe the turbulence observational data employed in this study. The data were retrieved from the Pacific Ocean (Fig. 1)[21]. They comprise profiles, each of which typically extends over a depth of 2000–6000 m below the sea surface, in turn comprising observational bins with width of . The turbulent energy dissipation rate for each bin, , is derived by averaging the observational values in the bin, which are estimated from the observed spectrum of the temperature vertical gradient based on the procedure presented in Goto et al.[21,22] (see Supplementary Information A for the estimation procedure). We restrict our investigation to the intermittency occurring at larger scales, .

Figure 1

Horizontal locations of the observed profiles (red) and land–sea boundaries (green). The units of longitude and latitude are and , respectively.

Let be the bin width, the horizontal coordinate of the i-th profile, and the vertical coordinate of the j-th point in the i-th profile. These positive-valued data exhibit the following characteristics:

Each profile defines an ordered set, which exhibits an extremely irregular evolution that impedes the recognition of a continuous curve along the depth direction (Fig. 2(a)).

Figure 2

Appearances of observed profiles.

After taking the logarithm of the values, the sequences appear to be more continuous (Fig. 2(b)).

If we normalise each value with the arithmetic mean along the profile to which it belongs as follows: then the histogram of the logarithmic values, appears as an asymmetric distribution, as we will see later in the Results section. Note the distinction between the two symbols; for the original energy dissipation rates, and for the normalised ones.

Multifractal analysis

We conduct a scaling analysis of the moments to derive the moment scaling exponent within the universal multifractal framework. Although the analysis could be extended to multidimensional objects[23], the limited number of samples (409 profiles) prevents us from conducting an extensive analysis in a multidimensional framework. Therefore, we treat each profile as an independent sample and analyse the statistical structure of the 1-dimensional object.

Universal multifractal model

The basic formulation of the universal multifractal model is as follows[23,24]: Suppose we have a multifractal field, , at resolution (), where r is the observational scale and L is the outer scale. The field is normalised by the mean, that is , which is conserved at all scales.

The probability of exceeding a scale-dependent threshold, , varies according to singularity aswhere represents the codimension function and represents equality up to multiplication by a slowly varying function of . Thus, the multifractal model is characterised by the property that the codimension varies with the singularity. This relation is equivalently represented as the scaling of the statistical moment of any order, q,where K(q) is the moment scaling function. The two functions, K(q) and , are actually related by the Legendre transformation because the moment generation function can be written in terms of the occurrence probability of singular events using the saddle-point approximation, [25], where . Functions K(q) and determine the variability of the multifractal field across the scales, .

Owing to a generalisation of the central limit theorem, several multiplicative processes comprising different generators converge to a universal multifractal[18,26], the moment exponent of which is expressed as follows:where is the multifractal index and is the codimension of the mean. Note that the case corresponds to the log-normal model advocated by the Russian school (Kolmogorov, Obukhov, Yaglom, etc.). This equation satisfies probability normalisation, , and energy conservation, . Its Legendre transformation gives:where .

Estimations based on the cascade model

In this section, we discuss the estimation of the energy input rate, , by utilising the information obtained from an observational profile. While the sample mean of the energy dissipation rate along a profile is simply indicated by the arithmetic mean of the vertical data values, the information on the energy input rate and its uncertainty still needs to be obtained. Therefore, we will estimate the posterior distribution of the energy input rate from observations. In particular, we focus on the median and confidence interval (CI). Although the arithmetic mean over a profile is the primary measure for the sample, the characteristics of the population can also be evaluated by using the joint probability density of several different sample statistics, obtained from the Monte Carlo simulation of the cascade model. The notation used in this section is summarised in Table 1.

Table 1

Notation for the estimation study.

Name	Notation	Definition
Index for vertical position	j	1,2,3,…,2^n
Energy dissipation rate	ϵ_j
Logarithm of energy dissipation rate	γ_j	logϵ_j
Energy input rate (or population mean)	ϵ¯
Logarithm of energy input rate	γ¯	logϵ¯
Median of estimated γ¯	γ_(0.5)	Pr(γ¯<γ_(0.5))=0.5
Stable Lévy generators	Γ_ik	∼S_α(σh^1/α,-1,-σ_α^^αh)
Width of Lévy generator	–	σh^1/α
Shift of Lévy generator	–	-σ_α^^αh = -σ^αcosπ2(2-α)h=-C_1α-1h
Logarithm of arithmetic mean	γ^	log2^-n∑j=12^ne^γ_j
Logarithm of geometric mean	γ~	2^-n∑j=12^nγ_j
Logarithm of quadratic mean	γ^♯	2^-1log2^-n∑j=12^ne^2γ_j
Marginal probability density function	q_1	Probability density of γ¯-γ^;
		q_1(·)=∫∫q_3(·,u,v)dudv
Joint probability density function	q_3	Probability density of (γ¯-γ^,γ^-γ~,γ^-γ^♯)

Multiplicative cascade simulation

To examine the relationship between various statistical quantities derived from observational profiles, we construct a simulation model for the multiplicative cascade by following the procedure described in Schmitt[27], as shown in Fig. 3. Each building block, , is a generator that obeys a left-skewed stable distribution, , with [28].

Figure 3

Schematic of the multiplicative cascade model. The energy dissipation rate at at resolution is considered as an example.

Consider a fixed horizontal position . Let be the energy input rate for a profile at , be the number of steps, be the scale index, and be the spatial index. The cascade simulation is performed for variable as follows. The output, , represents the logarithm of the energy dissipation rate at the horizontal position, , and the vertical position, , at the resolution . By using the floor function, the cascade process can be more compactly represented asAn important implication of this formulation is that the arithmetic mean of the vertical data points is not necessarily equal to the the energy input rate because the cascade process has a fluctuating nature. In other words, a realisation of the vertical average, , is not always equal to , whereas the expectation is; hence, we can regard the latter as the population mean for a profile. Below, we focus mainly on the relationship between the arithmetic mean over a profile and the energy input rate. We perform statistical estimations from one to the other of these quantities based on the cascade model.

For each spatial index , set .

For each scale index , repeat the following steps:

For each spatial block , perform the following steps:

(a)  Generate a random variable, , which obeys [29].

(b)  For each spatial index , downscale X by

For each spatial index , set .

Estimation of the energy input rate

In this subsection, we first describe the statistical relationship between the population mean and various moments over a profile, based on the cascade model. Then, we derive a formula for the posterior probability given the observation of the moments. Finally, we use this formula as the basis of a concrete procedure for computing the posterior probability.

We assume that each set of ’s is generated by an n-step cascade model as in Eq. (12). Here, we want to estimate the energy input rate , which corresponds to the population mean, by using the information from the observed data . In this regard, in addition to the arithmetic mean , which corresponds to K(1) in Fig. 5, we can also use other moments over a profile, e.g., the geometric mean and quadratic mean , which correspond to and K(2), respectively. We can derive the following expressions based on Eq. (12):where we find that the term is factored out. Therefore, , , and are independent of , and thus dimensionless.

Figure 5

Moment scaling exponent K(q) for observational data (cyan). Best-fitting multifractal model with stable Lévy generators (black), and with Gaussian generators (red). Each error bar in cyan shows the standard deviation for the fitting of K(q). Dotted lines in black and red indicate the ranges of error due to the uncertainty of parameters in the corresponding models.

The structure of the cascade model implies that the appearance probability of given is determined only by their difference: . Furthermore, by assuming that we have no prior information about , Bayes’ theorem is applied to invert it into the posterior probability for as follows.Furthermore, we can extract information from and . They are encoded in the joint probability density function (PDF) computed from Monte Carlo simulation of the cascade model. Then, a conditional PDF is derived asA procedure similar to Eq. (16) can be applied to obtain another posterior probability:under the constraints .

On the basis of the above formulation, we perform an identical twin experiment obeying the following procedure. The same procedure is applied to the real data experiment, except that in 2 is unknown, as follows. Among the indices for the estimated result, the conditional expectation of is not necessarily defined as a finite value because the posterior distribution of is neither Gaussian nor left-skewed stable. In contrast, the percentiles, including the median and the CI, are always defined for the distribution. They are also preserved, , under the increasing transformation . We therefore employ the median and the CI as robust indices.

Perform a Monte Carlo experiment to obtain .

Set the energy input rate to .

Create many random samples of the profile using the cascade model in Eq. (12).

Add up the frequency of occurrence to derive the joint PDF .

Set the energy input rate to a random number.

Create a random pseudo-observation sample of the profile using the cascade model in Eq. (12).

Calculate the statistics for the profile.

Compute the conditional PDF .

Calculate the median and CI for the estimated .

Compare the estimate of with its true value.

Pick an observed profile, and calculate the statistics for the profile.

Compute the conditional PDF , using the joint PDF obtained from the Monte Carlo experiment.

Calculate the median and CI for the estimated .

To confirm that the uncertainty in the estimated values of and does not diminish the performance of the proposed method, we treat these parameters as random variables (specified in the results), when making pseudo-observation samples in the identical twin experiments (procedure 3).

We have thus established a procedure for estimating the population mean for a profile based on the joint probability distribution of several moments over a profile computed by a Monte Carlo simulation of the cascade model.

Results

Analysis of observational data

Suppose we have the observational data of the normalised energy dissipation rate, , in bin width at the horizontal position , as well as their spatial average in width . In terms of the universal multifractal model (8), the scaling of the statistical moments in the data takes the formwhere denotes the expected value. This implies that the expectation of the q-th moment at a scale over the one at another scale should be equal to the K(q)-th power of the resolution ratio, regardless of the horizontal position . We approximate the expected value in Eq. (19) with the empirical averagewhere is a resolution larger than or equal to , and is the normalised value defined in Eq. (5). The superscript (i, k) runs across all profiles indexed by i, each of which has total segments in resolution r. By substituting Eq. (20) into Eq. (19), we can evaluate the values of K(q) according to q.

The scalings for several moments are shown in Fig. 4. Using the various slope values, the observational curve of (q, K(q)) in the range of is indicated in Fig. 5 in cyan.

Figure 4

Scale dependency of the moments, where is the width of the observational bin. The moment scaling exponents are found to be .

We can estimate the parameters, and , by fitting the theoretical curve (9) to the observational curve. To consider the uncertainty in the observational curve, we used the bootstrap method[30] with 1000 trials, each of which has 409 profiles that are randomly sampled with replacements from the original set of 409 profiles. We thereby obtained the parameters for the multifractal model with stable Lévy generators, i.e., the universal multifractal model. By taking into account the dependency on , we can also estimate as . By a similar procedure, we obtain the parameter for the multifractal model with Gaussian generators, corresponding to the original Yaglom cascade with . In Fig. 5, the observational curve (cyan) and the theoretical curve for the multifractal model with stable Lévy generators (black) are in good agreement, while the theoretical curve for the multifractal model with Gaussian generators, i.e., the log-normal model (red), has a different curvature from the observational curve.

The parameter values for the multifractal model with stable Lévy generators are largely consistent with previous results for atmospheric dissipation fields ( for the horizontal shear of a velocity field[31]; for vertical kinetic energy flux[32]).

Figure  6 shows the theoretical curve of extremes for the multifractal model (10) in black and the observational curve,in cyan, where is used; this is a typical scale ratio in the data. Note that the correction term,compensates for the prefactor in the asymptotic complementary cumulative distribution function, [28], Eq. 1.2.11]. The two curves (10) and (21) appear to be in general agreement, except for a slight discrepancy that is possibly due to the ambiguity in the selected scale ratio, . Moreover, as our data have the sampling dimension[23] , the upper bound for q is calculated to be (the slope of the navy-blue line in Fig. 6), which justifies the range we set ().

Figure 6

Codimension of singularities for the best-fitting multifractal model with stable Lévy generators (black). The corresponding curve for the observational data is shown for reference (cyan). Sampling dimension and the limitation for the moment exponent (the slope of the navy-blue line) are also shown. Dotted lines in black and navy indicate the ranges of error due to the uncertainty in the model parameters.

To demonstrate the appropriateness of the universal multifractal model, the histogram for the logarithm of the bin values in the observational data is shown in Fig. 7 and compared with the samples from multiplicative cascade models. Each bin value is normalised by the arithmetic mean along the profile it belongs to: The histogram for the logarithm of bin data, , appears to be in good agreement with the histogram of samples generated by the 8-step cascade model with stable Lévy generators (black; ) and in poor agreement with that generated by the multifractal model with Gaussian generators, i.e., the log-normal model (red; ).

Figure 7

Distribution of the logarithm of observational data normalised for each profile (cyan), and comparison with the statistics of samples generated from multiplicative cascade with Gaussian/stable Lévy generators (red/black).

Moreover, in the same manner as the correlation in Yaglom’s cascade, the observational profiles have a power-law autocorrelation,where is the energy dissipation rate at depth z. The negative exponent explains the discontinuous characteristics observed in the profiles (see Fig. 2a).

Simulations of cascade model

Before estimating the energy input rate that corresponds to each observational profile, we first performed a Monte Carlo experiment with simulations of particles (profiles) using the 8-step cascade model. The constants and parameters used in the simulation and estimation study are summarised in Table 2.

Table 2

Constants and parameters for the estimation study.

Meaning	Parameter	Value
Number of steps in cascade	n	8
Number of vertical points	2^n	256
Step size of cascade	h	log2
Number of samples in Monte Carlo simulation	M	1.024×10^10
Multifractal index	α	1.62
Codimension of the mean	C_1	0.352
Width of bins in histogram of q_3	(w_1,w_2,w_3)	(0.1, 0.1, 0.05)
Number of bins in histogram of q_3	(n_1,n_2,n_3)	(500, 200, 200)
Number of profiles in identical twin exp.	M^′	30,000
Number of observed profiles (Obs.)	I	409
Obs. with more than 2^n vertical points	–	353

For each particle (or profile), we generate random numbers from according to the procedure discussed in Methods, and we add up the histograms for all the particles into the joint PDF (Fig. 8).

Figure 8

Examples of cross section cut of joint probability density function . (a) section cut with section lines , and (b) section cut with section lines .

Examples for conditional PDF are shown in Fig. 9.

Figure 9

Examples of conditional probability density function (coloured) and marginal probability distribution (black), along with median (dashed lines) and confidence intervals (dotted lines). The conditional probabilities are for .

Identical twin experiment

Second, we estimated the posterior probability distribution for the energy input rate, by inverting the probability distribution computed from the simulation of the cascade model. Before applying this to real data, we performed an identical twin experiment using pseudo-observational data, whose energy input rates were given manually; thus, the estimation result could be checked against them. The inversion was performed using the result of the joint PDF created by the cascade model with the Lévy generator (case L3). When generating the pseudo-observational data, the parameters for each profile were set randomly as using standard normal random numbers, and .

The result of the identical twin experiment using is shown in Fig. 10a. In 28,497 trials out of 30,000 (about ), the true value of lies within the CI, which ensures the validity of the estimation method.

Figure 10

Result of identical twin experiments. Median (horizontal axis) versus confidence interval (vertical axis) of . (a) Result using probability density function . (b) Result using probability density function . For given values of the median on the horizontal axis, the points on the vertical axis indicate the values of the confidence interval (purple segment), arithmetic mean (orange), geometric mean (blue), and energy input rate (black). For readability, 300 points and 60 intervals are drawn out of 30,000 trials.

Control experiments

To show the superiority of the proposed method (L3), which performs an inversion using the result of the joint PDF created by the cascade model with the Lévy generator, the generated samples above are also estimated by other methods: an inversion based on the marginal PDF , an inversion based on the joint PDF generated by a multiplicative cascade with Gaussian generators, and a simple bootstrap method.

Error estimation based on the PDF of the arithmetic mean

Using several different moments for a profile should be effective for precise error estimation. To check this, we perform a control experiment using the marginal PDF, , of the arithmetic mean (case L1). We take samples created by the cascade model with the stable Lévy generators, and then estimate the CI via the marginal PDF . The results are shown in Fig. 10b, where the CI and median show a common positional relation to the arithmetic mean. Comparing the case with to the one with , of trials in the former have narrower CIs. This indicates that using information from and improves the error estimation.

Error estimation based on a cascade model with Gaussian generators

Stable Lévy generators have asymmetry in the distribution, which also affects the mean and median. Therefore, it is necessary to use such asymmetric generators for the error evaluation. To check this, we also evaluate the samples in the twin experiment using the probability distribution generated by a cascade model with Gaussian generators for comparison. We take samples created by the cascade model with stable Lévy generators (case G3) or Gaussian generators (case G1), and then estimate the CI via the joint PDF based on the cascade model with the best-fitted Gaussian generators. The results are shown in Fig. 11a, where a significant portion (9891 trials out of 30,000) of the true energy input rate (black) protrudes above the CI (purple). This means that the values of the energy input rate are underestimated if we assume a Gaussian distribution for the generators, and it also illustrates that it is inappropriate to use the statistics from the simulations of a cascade model with Gaussian generators for error evaluation.

Figure 11

Result of control experiments. Median (horizontal axis) versus confidence interval (vertical axis) of , which are obtained using (a) joint probability density function based on the cascade model with Gaussian generators, and (b) the bootstrap method. For given values of the median on the horizontal axis, the points on the vertical axis indicate the values of the confidence interval (purple segment), arithmetic mean (orange), geometric mean (blue), and energy input rate (black). For readability, 300 points and 60 intervals are drawn out of 30,000 trials.

Error estimation using the bootstrap method

The simplest method for error evaluation is applying the bootstrap method to each profile. However, the errors cannot be properly assessed by such a conventional method. To verify this, we estimate the CIs by applying the bootstrapping method to the twin experiment. For each trial, we use a 1000-member ensemble for the bootstrapping. Each member is constructed as follows: If a profile at the horizontal point has a set of observations of the energy dissipation rate, , we randomly take samples with replacement from the set. The results of the error evaluation of the mean dissipation rate are shown in Fig. 11b. The CIs are evaluated very narrowly, and in many trials (22,164 trials out of 30,000), the true energy input rate (black) is outside the CI, which indicates that the error estimate is far too optimistic. This illustrates that it is irrelevant to use the conventional bootstrap method for the error evaluation.

Comparison of skill

For a fair comparison of skill in the above control experiments, we examine the errors for an estimator of the energy input rate. First, we define the estimator, , from the median estimate, , aswhere [m] represents the m-th sample, and is the median of the samples generated from the energy input rate . Here, a is defined to empirically satisfy the unbiasedness:Then, in terms of , we define the relative error with sample size asNote that the error might not necessarily converge when because the distribution of is not Gaussian. Nevertheless, we can still evaluate the error for a finite and use it for the comparison of skill.

The errors with for various conditions are listed in Table 3. The combinations of the cascade model with the Lévy generator or Gaussian generator and the use of joint PDF or marginal PDF are compared; these correspond to cases L3, L1, G3, and G1 above. The error for the arithmetic mean is also shown for reference. The smallest error among these is for the estimator using the result of joint PDF created by the cascade model with a Lévy generator (L3). The errors for cases L1 and G1 are comparable to that for the arithmetic mean because these estimators are constructed from the statistics of the arithmetic mean. This result clearly shows that the proposed method (L3) can be used to define an estimator that yields superior estimates of the energy input rate than the arithmetic mean or other methods (L1, G3, or G1).

Table 3

Comparison of the error, , for various estimators in identical twin experiment.

ind_M^′	q_3		q_1
Lévy	L3	0.89 (0.745)	L1	0.98 (0.718)
Gauss	G3	1.0 (1.61)	G1	0.98 (0.689)
γ^	0.99 (1.00)

Lévy/Gaussian indicates the generator used in cascade model simulation. Here, / indicates whether joint PDF or marginal PDF is used as the density generated by the cascade model simulation; indicates the arithmetic mean. Each number in parentheses is the prefactor a for the corresponding estimator.

Real data experiment

We applied the same procedure as in the identical twin experiment to the real observational profiles of the energy dissipation rate. Each profile was characterised by , , and , which were utilised as observational constraints. By means of inversion, we derived the CI of for each profile at different horizontal locations. The estimated CIs for real data are shown in Fig. 12b; these CIs exhibit a similar appearance to the ones for the identical twin experiment in Fig. 12a. Figure 13 shows the estimate for CIs on the sections along and . Along , the median of rarely exceeds , except around , , or . The peak of the arithmetic mean at is approximately 2.5 times larger than the median estimate, which can lead to overestimation.

Figure 12

Results of the real data experiment compared with those of the identical twin experiment. Median (horizontal axis) versus confidence interval (vertical axis) of . (a) Results of identical twin experiment. (b) Results of real data experiment. For given values of the median on the horizontal axis, the points on the vertical axis indicate the values of the confidence interval (purple segment), arithmetic mean (orange), geometric mean (blue), and energy input rate (black). In b, only 70 confidence intervals out of 353 trials are shown for readability.

Figure 13

Geographical distribution of median (purple dots) and confidence interval (purple shade) of along (a) and (b) . The horizontal axis shows the location. Arithmetic mean (orange) and geometric mean (blue) are also shown.

Along , the median shows several significant peaks over at around , , and 27 to . We could have underestimated the peaks at around and 27–29, but overestimated the one at around , if only the arithmetic means were used.

For the analysis of the observations on the section along , we should consider the effects of repeated observation. In fact, the observations were performed twice at some horizontal locations. For such cases, we simply assume that two independent realisations of a common are observed. In this regard, the inversion formula in Eq. (16) is modified as follows:where two observations are distinguished by the subscripts 1, 2. This distribution is the normalised product of the two distributions. When considering and , the same procedure as in Eq. (28) is applied to instead of :Figure 14b shows the estimation of CI along by taking into account the effect of repeated observation. For comparison, the result without considering the repeated observation is shown in Fig. 14a, where each observation is assumed to correspond to independent . We can see that the CIs become narrower when considering the effect of repeated observation. Furthermore, along , the median shows a significant plateau on the order of at around , and a significant peak on the order of at around .

Figure 14

Geographical distribution of median (purple dots) and confidence interval (purple shade) of along . The effect of repeated observation is considered in (b), but not in (a). The horizontal axis shows the location. Arithmetic mean (orange) and geometric mean (blue) are also shown.

Discussion

We have analysed the observed data obtained from oceanic turbulence measurements and shown that the vertical sequences of energy dissipation rates in a profile have an intermittent structure that obeys a scaling law. In this study, we have laid greater emphasis on the ’energy input rate’, which refers to the population mean for a profile, than on the ’mean energy dissipation rate’, which is the sample (arithmetic) mean over a profile. Based on the scaling property, we have proposed a method of estimating the energy input rate, given the sample statistics of an observed profile.

For scaling within the observed profiles, the statistical properties of our data are consistent with the universal multifractal model, which has a moment scaling exponent of with a multifractal index and codimension of the mean . This result elucidates the universality that is inherent in the vertical structure of oceanic turbulence data.

The energy input rate and its uncertainty can be estimated using the results of Monte Carlo simulation of the cascade model with stable Lévy generators. This method computes the conditional probability, given the observed values of the arithmetic mean, geometric mean, and quadratic mean over a profile. The estimate provides additional information on the uncertainty of the energy input rate.

Furthermore, a comparison to control experiments has demonstrated that the proposed method is superior to a simple bootstrap method, an inversion based on the PDF generated by a multiplicative cascade with Gaussian generators, or an inversion based on the PDF of the arithmetic mean.

A real data experiment using the observed profiles has demonstrated a geographical distribution of the median estimates and confidence intervals of energy input rate, providing information on the range of values in which the turbulent energy can be dissipated per unit depth at each horizontal location.

Thus, we have found an answer to the question: ‘How can one estimate energy input rate from the vertical profile data of the energy dissipation rate?’ By analysing the intermittency in the observed data, we can construct a multiplicative cascade model based on the universal multifractal formalism that can reproduce the statistics of the data. Then, based on the observed data, the energy input rate can be estimated by inverting the probability distribution obtained from Monte Carlo simulations of the cascade model.

Since the observed sequence of the energy dissipation rate fluctuates greatly due to intermittency, it is difficult to extract robust information from an observed profile only by examining the average. Using the proposed method, it is possible to estimate the population mean for a profile even when repeated observations cannot be made at the same horizontal position. In other words, we can distinguish, to some extent, whether a profile shows an occasional large mean or whether the population mean itself is large. Therefore, more information can be extracted from a small amount of turbulence observation data, which can be a great advantage in regional data analysis.

Theoretically, this technique can easily be extended by utilising more statistics over a profile besides the arithmetic mean, geometric mean, or quadratic mean. Note however that such an extension may easily suffer from the curse of dimensionality, and thus, it can become impractical.

Even though we have used a discrete cascade model for simplicity and computational viability, we can extend it to a continuous cascade[33], which may improve the estimation accuracy at the cost of increased computational burden.

To investigate the scaling of the velocity spectrum in the horizontal and vertical directions, we should apply anisotropic scaling theory (e.g., the Kolmogorov–Bolgiano–Obukhov model[23]). However, this study only aimed at investigating the scaling of the energy dissipation rate in the vertical direction in a purely statistical manner. Nevertheless, because the energy dissipation rate is one of the key quantities in scaling analysis, our results will contribute to further studies on how intermittency affects various scaling behaviours of turbulence in buoyancy-driven stratified fluids.

Supplementary Information.