A variable Prandtl number power spectrum of optical turbulence and its application to oceanic light propagation

Light interaction with the turbulent ocean can be fully characterized with the help of the power spectrum of the water's refractive index fluctuations, resulting from the combined effect of two scalars, temperature and salinity concentration advected by the velocity field. The Nikishovs' model [Fluid Mech. Res. 27, 8298 (2000)] frequently used in the analysis of light evolution through the turbulent ocean channels is the linear combination of the temperature spectrum, the salinity spectrum and their co-spectrum, each being described by an approximate expression developed by Hill [J. Fluid Mech. 88, 541562 (1978)] in the first of his four suggested models. The forth of the Hill's models provides much more precise power spectrum than the first one expressed via a non-linear differential equation that does not have a closed-form solution. We develop an accurate analytic approximation to the forth Hill's model valid for Prandtl numbers in the interval $[1,1000]$ and use it for the development of a more precise oceanic power spectrum. To illustrate the advantage of our model we include numerical examples relating to the spherical wave scintillation index evolving in the underwater turbulent channels with different average temperatures, and, hence, different Prandtl numbers for temperature and salinity. Since our model is valid for a large range of Prandtl numbers it can be readily adjusted to oceanic waters with seasonal or extreme average temperature and/or salinity or any other turbulent fluid with one or several advected quantities.


Introduction
The growing interest in the development of optical underwater imaging systems [1], laser communications [2,3] and remote sensing schemes [4] stipulates the request for the occurate predictions of the oceanic turbulence effects on the light wave evolution.Such effects can be characterized with the help of the parabolic equation, the Rytov approximation, the extended Huygens-Fresnel integral and other classic methods as long as the refractive-index power spectrum of turbulent fluctuations is established [5,6,7,8].As applied to oceanic propagation, the theoretical predictions for the major light statistics of the optical waves have been established (see recent review [9]) on the basis of the widely known power spectrum developed by Nikishovs [10].In particular, the spectral density and the beam spread for coherent and random beams have been analyzed [11], the spectral shifts in random beams were revealed [12], the polarimetric changes of the random electromagnetic beams have been discussed [13], the loss of coherence of initially coherent beams has been addressed [14], the scintillation analysis was provided in Refs.[15] and [16] and the structure functions of various waves were derived in Refs.[17] and [18].
In the situations when there is more than one scalar field advected by turbulence, the total refractive-index power spectrum can be expressed as a linear combination of the individual spectra and their co-spectra.For instance, in the case of the oceanic turbulence with advected temperature and salinity (sodium chloride) concentration, the total spectrum is the linear combination of three spectra: temperature spectrum, salinity spectrum, and their co-spectrum [10] (see also [19]).An important dimensionless quantity determining the turbulent spectrum profile is the Prandtl number, P r, defined as a ratio of viscous diffusion rate (kinematic viscosity) [m 2 s] to the thermal diffusion rate (thermal diffusivity) [m 2 s].Since the Prandtl numbers for temperature and salinity in the water are sufficiently large (greater than one) their power spectra have three regimes: inertial-convective regime with the Kolmogorov power law Φ n (κ) ∼ κ −11/3 at the low wave numbers, viscous-convective regime with the power law Φ n (κ) ∼ κ −3 for higher wave numbers and the viscous-diffusive regime at even higher wave-numbers where the spectra decrease very rapidly driven by the diffucsion of the advected quantity.Due to two different power laws at the high wave numbers each spectrum forms a characteristic "bump" [20,21,22].
To achieve the precise match of the spectra at the bondaries between different regimes one of the four models proposed by Hill [23] may be applied.Among them, the Hill's model 1 (H1) that provides the simple analytic approximation but lacks prescision in accounting for the bump, and the Hill's model 4 (H4) that is precise but involves a non-linear differential equation without known analytic solution.‡ The widely used Nikishovs' oceanic turbulence spectrum [10] is developed as a linearized polynomial involving temperature spectrum, salinity spectrum and their cospectrum, each being based on H1 see (Appendix I for more details).Hence, it has a simple analytic form but lacks precision.It can, in principle, be applied to various Prandtl numbers but, as we will illustrate, substantially deviates from H4, especially for Prandtl numbers much larger than one.Also, a recent model for the oceanic power spectrum [22] was obtained by fitting H4 at specific P r values: P r = 7 for temperature spectrum, P r = 700 for salinity spectrum and P r = 13.86 for their co-spectrum.It was ‡ Further, H4 agrees with the Kraichnan spectrum at high wave numbers, while H1 evidently deviates from it.The precision of the Kraichnan spectrum and, hence, of H4 is well established [24,25,26,32].revealed that the oceanic power spectrum of [22] has a substantially higher bump at high spatial frequencies as compared with that for the Nikishov's spectrum, resulting in a greater scintillation index in light waves interacting with such a medium.
Moreover, within the last two decades, several mathematically tractable models based on H4 have also been introduced for atmospheric refractive-index spectrum [28,29,30,31,27], by fitting P r = 0.72 and P r = 0.68 for temperature and humidity Prandtl numbers, respectively.For such low Prandtl numbers in the atmospheric turbulence the viscous-convective regime is not that prevalent but nevertheless the highfrequency bump can still be formed.Among the precise atmopsheric power spectrum models are the Frehlich's model being a fifth-order Laguerre expansion [29], as well as the Grayshan's and the Andrews' models both being products of polynomials and Gaussian functions [31,30].All these power spectra provide convenience for theoretical calculations, suggesting that the precise analytical fits do not have to be mathematically involved.
The aim of this paper is to employ the H4 model for developing a power spectrum that can be readily adapted to the wide range of Prandtl numbers but, at the same rate, remain numerically accurate and mathematically tractable.Such new spectrum can be applied for classic homogeneous oceanic turbulence in a variety of thermodynamic states.Moreover, it can also be employed for optical characterization of other turbulent media, with one or two advected quantities, for example, fresh turbulent water contaminated by a chemical.The ratio of kinematic viscosity to diffusivity is orders of magnitude larger for temperature and salinity fluctuations in water than for temperature and humidity fluctuations in air [21], resulting in much larger Prandtl numbers (averaging at P r T = 7 for temperature and P r S = 700 for salinity [9,10,20,22].More importantly, these numbers can exhibit variation.For example, at 1 Bar pressure (water surface), at fixed salt concentration at 35ppt, as the temperature increases from 0 o C to 30 o C, P r T decreases from 13.349 to 5.596 and P r S decreases from 2393.2 to 456.1 (more details in Appendix II).On the other hand certain oil-like substances can lead to Prandtl numbers up to 10 5 .
The paper is organized as follows.In section 2.1 we introduce an approximate fit to H4 for P r ∈ [1, 1000] for handling situations in which a single scalar is advected by turbulence.Section 2.2 provides the validity estimation of the fit by means of the dissipation constraint.Section 3.1 introduces an H4-based power spectrum of oceanic turbulence by applying the fit of Section 2.1 to three spectra: temperature-based, salinity-based and the their co-spectrum.In Section 3.2 we apply the new spectrum model for analytic calculation of the scintillation index of the spherical wave for a variety of Prandtl numbers and illustrate its agreement with numerical calculation based on the H4 model as well as its descrepancy with the result based on the Nikishovs' spectrum.In Section 4 all the results are briefly summarized.Let us begin by developing the three-dimensional (3D) scalar power spectrum of the refractive index fluctuations produced by a single quantity advected by a fluid, as an analytical fit to the H4 model.According to [32] the 3D scalar spectrum Φ n (κ) is related to a universal dimensionless function g(κη) by expression where β is the Obukhov-Corrsin constant; ε is the dissipation rate of the kinetic energy [m 2 s −3 ]; η is the Kolmogorov microscale [m]; χ is the ensemble-averaged variance dissipation rate.The units of χ depend on the advected quantity, for instance, for temperature, salinity and refractive index they are K 2 s −1 , g 2 s −1 and s −1 , respectively.According to model H4, g(κη) is the solution of the nonlinear ordinary differential equation [23,32] d dz (z 2b + 1) where and z = κηa −1 , c = βa 4/3 P r −1 .Generally, as recommended by Hill [20] the following fixed values of parameters a = 0.072, b = 1.9, and β = 0.72 § will be used thoughout the paper.The numerical solution g(κη) of Eq. ( 2) exhibits an obvious "bump" with height g max gradually increasing with growing values of P r, for P r > 0.2 [32].Considering the significant effect of P r on g(κη) [32], we have set the form of g(κη) as where h 2 , h 5 and h 7 are limited to positive numbers and then fitted the numerical solution g(κη) of Eq. (2) with Eq. ( 4) by iterative least squares method.This resulted in the following expression for g(κη): g(κη) = 1 + 21.52(κη) 0.63 c 0.02 − 18.17(κη) 0.58 c 0.04 × exp −176.41(κη) 2 c 0.96 .
(5) § The value of β may change with environment, and has always been used as an experimental constant.Several experiments gave its range 0.22 ∼ 0.78 [37].The value 0.72 has been widely used in oceanic optical research, so we use the value 0.72 in Table 1 and Figure 2, but keep it as a parameter in all formulae.
valid in the interval P r ∈ [1, 1000].On substituting from Eq. ( 5) into Eq.( 1), we find that the fitted power spectrum based on H4 takes form Φ n (κ) = 1 + 21.52(κη) 0.63 c 0.02 − 18.17(κη) 0.58 c 0.04 To examine the prescision of fitted solution g(κη) we compare it with analytic formula H1 (see Appendix I) and H4 (calculated numerically) for several values of the Prandtl number.Figure 1(a) produced for P r = 1 implies that the fitted function g(κη) calculated from Eq. ( 5) and that obtained from H1 model agree with that based on H4 model reasonably well while Eq. ( 5) provides with a better bump-shape reconstruction.Further, figures.1(b)-(d), obtained for P r = 10, 100 and 1000, respectively, illustrate that g(κη) in Eq. ( 5) is very consistent with that based on H4, but it is not the case for g(κη) based on H1 the latter leading to drastic underestimation of the bump's strength.Hence, if one uses H4 as a benchmark, the proposed g(κη) in Eq. ( 5) and, hence, the Φ(κ) in Eq. ( 6) show advantages in describing the spectral bump for all P r ∈ [1, 1000], and especially so for P r 10.
We note that g(κη) in Eq. ( 5) has the same functional form as those in the Andrews' [30] and in the Grayshan's [31] power spectra being the product of a polynomial and a Gaussian function.On the other hand, it also relies on κ and P r independently, making the bump height g max increase with growing P r values, which is of importance in achieving the prescision in the spectrum model valid for the wide range of P r.

Dissipation constraint analysis
Although a good numerical fit for the power spectrum was obtained in the previous section for a wide range of P r it appears of importance to estimate whether the new model agrees with the basic laws of fluid mechanics.In this section we will verify to which degree the dissipation constraint [29,27] being the consequence of the scalar transport equation is satisfied by Eq. ( 5).The dissipation constraint can be expressed via integral [29]: The closer the value of X to 1 is more g(x) agrees with the dissipation constraint.
The values of X calculated on substituting Eq. ( 5) into Eq.( 7) are shown in Table 1 for different values of the Prandtl number.The table entries imply that Eq. ( 5) underpredicts X by 8.0% and 17.3% when P r = 1 and 1000, respectively.This descrepancy is intensified with the increasing Prandtl number which, in its turn, corresponds to the increasing bump height.This indicates that a more precise fitting of the spectral bump occurs at the expense of deviation from the dissipation constraint.For comparision, X corresponding to other models were calculated in Ref. [27], and are shown in Table 2.These models fit H4 at P r = 0.72.Apart from the Frehlich model the values of X for the rest of the models do exhibit certain deviation from value 1.On comparing entries of Tables 1 and 2 it appears that the model in Eq. ( 5) meets the dissipation constraint in a satisfactory manner.

Models
Frehlich Churnside Andrews Grayshan X for P r = 0.72 1.00 1.12 0.95 1.37 The power spectrum model given by Eq. ( 6) is the first of the two main results of our paper.It can be used for any turbulent medium in which a single scalar is advected by the velocity field and is applicable for the wide range of the Prandtl numbers, P r ∈ [1,1000].Moreover, it is confirmed that model ( 6) is in the reasonable agreement with the dissipation constraint.

Practical example: oceanic refractive-index spectrum
In this section we will demonstrate how to extend the power spectrum model given by Eq. ( 6) from one to two advected scalars with different Prandtl numbers.The best illustration of such a medium is the oceanic turbulence, hence we will adopt the corresponding notations to maintain practicality.
On following the Nikishov's linearized polynomial approach [10] we assume that the fluctuating portion of the refractive index n is given by linear combination of the temperature fluctuation T = T − T and the salinity concentration fluctuation S = S − S , T and S representing the instantaneous temperature and salinity values and the angular brackets standing for the statistical average.In Eq. ( 8) A is the thermal expansion coefficient and B is the saline contraction coefficient [10].Therefore, the power spectrum of the refractive index fluctations can be expressed as the following linear combination of the temperature spectrum Φ T (κ), the salinity spectrum Φ S (κ), and the co-spectrum Φ TS (κ) [10]: However, unlike in Ref. [10] where each of these three spectra was based on the H1 model, we will use the fit for the H4 model obtained above [see Eq. ( 6)], i.e., Φ i (κ) = 1 + 21.52(κη) 0.63 c i 0.02 − 18.17(κη) 0.58 c i 0.04 where c i = 0.072 4/3 βP r −1 i , P r T and P r S are the temperature and salinity Prandtl numbers, respectively, P r TS = 2P r T P r S (P r T + P r S ) −1 is the coupled Prandtl number [19,22], and χ i are the ensemble-averaged variance dissipation rates which are related by expressions [33] d r being the eddy diffusivity ratio and ω being the relative strength of temperaturesalinity fluctuations.In most practical cases, d r and ω are related as [33] The power spectrum model given by Eqs. ( 9) and ( 10) is the second of the two main results of our paper.It gives the precise analytic description of optical turbulence with two advected quantities while each of them may have P r ∈ [1,1000].In its present form the new model is specifically applied for the oceanic waters in which the turbulent velocity field advects temperature and salinity concentration.However, it can be readily adapted for a variety of other turbulent media based on two scalars, for instance a fresh water contaminated by a chemical substance or any exotic turbulent fluid.Futhermore, a straightforward linearization of three or more advected quantities can be worked out in a similar manner.

The scintillation index of a spherical wave
We will now illustrate how the analytical fit for the oceanic power spectrum given by Eqs. ( 9) and ( 10) compares with the Nikishovs' spectrum (based on analytical model H1) and the H4 model (handled via numerical calculations) as it manifests itself in the evolution of the scintillation index of light propagating though an extended turbulent channel.The scintillation index of an optical wave at distance L from the source plane is generally defined by expression [5] where I is the instantaneous intensity of the wave, and the angular brackets stand for the statistical ensemble over the medium realizations.We will restrict ourselves to the scintillation index of a spherical wave being one of the most important parameters of light-turbulence interaction.In addition to the Rytov variance, i.e., the scintillation index of the plane wave, it can serve as an indicator of the global turbulence regime (weak/focusing/strong) [5].This quantity was extensively studied for the atmopsheric propagation in dependence from a number of turbulence parameters [5].Based on the Nikishovs' spectrum the scintillation index of the spherical wave was also previously evaluated on propagation in the oceanic water (c.f.[12]).
Let us now analytically derive the expression for the scintillation index of the spherical wave in the water channel described by spectrum in Eqs. ( 9) and (10).According to the Rytov perturbation theory it takes form with where i = √ −1 (not to be confused with subscript i also used in the paper), k = 2πn 0 /λ 0 , λ 0 is the wavelength in vacuum, n 0 (ca.1.33) is the averaged refractive-index of the ocean water.
On substituting from Eq. ( 9) into Eqs.( 14) and ( 15) we arrive at the scintillation index in of the form and possibly variable Prandtl numbers the accurate analytical fits have not yet been introduced.
In this paper we have introduced a very simple yet occurate analytical model for

2 .
Analytic Approximation of the Hill's Model 4 for One Advected Quantity 2.1.The Polynomial-Gaussian fit of the Hill's Model 4

Table 1 .
Values of X calculated from Eq. (5) for various Prandtl numbers.

Table 2 .
The values of X corresponding to reported models