This article examines some of the difficulties associated with the determination of C2n over water in a coastal region using a bulk model. The analysis shows the need to supplement bulk models with elements that do not belong to traditional Monin–Obukhov surface-layer theory. A reexamination of the scintillation measurements collected during the Electro-Optical Propagation Assessment in a Coastal Environment (EOPACE) campaign leads the authors to include 1) the nonuniformity of sensible heat and humidity fluxes and 2) a deficit of the scintillation that seems to depend only on the characteristic virtual potential temperature. It is suggested that the anomalous scintillation deficit represents an alteration of a characteristic length related to the optical turbulence, likely caused by the interaction of the surface layer with the sea surface. The new parameters are estimated using Bayesian regression methods applied to the EOPACE data. Predictions obtained with this new model are then compared with scintillation measurements obtained during the recent Validation Measurement for Propagation in the Infrared and Radar (VAMPIRA) campaign. Better agreement is obtained with the new model than with a conventional bulk model. The implications of the modifications made to the calculation of C2n are discussed.
Atmospheric turbulence can significantly affect the performance of electro-optical imaging systems to detect and identify objects, being responsible for scintillation, blurring, and image wandering. The theory, in the scientific literature, relates the turbulence effects to the optical refractive index structure parameter C2n, which characterizes the rapid fluctuation of the index of refraction (Beland 1993). For predicting turbulence effects on an imaging system, it would be beneficial to use valid bulk estimates of C2n obtained from standard environmental measurements.
Bulk models of C2n based on Monin–Obukhov similarity (MOS) theory were originally developed to describe turbulence over land (Andreas 1988). For calculations above the sea, new terms incorporating the effect of humidity, such as in the models developed by Claverie et al. (1995) and Kunz et al. (1996), have been proposed. These models were designed to take into account humidity effects made apparent by studies of optical scintillation and bulk meteorological measurements made over the sea (Friehe 1977; Davidson et al. 1981; Kunz et al. 1996). In their analysis of transmission measurements conducted during the Electro-Optical Propagation Assessment in a Coastal Environment (EOPACE) field experiment over the bay in San Diego Bay, California, Frederickson et al. (2000, hereinafter FDZB) showed for the first time an extensive comparison between scintillation measurements above water and bulk predictions. Their study revealed that bulk predictions can strongly underestimate the observed scintillation when the air–sea temperature difference (ASTD) is close to zero and overestimate them when the magnitude of ASTD is appreciable. More recently, scintillation measurements were performed in a different environment at the Validation Measurement for Propagation in the Infrared and Radar (VAMPIRA) trial conducted over water in Eckernförde Bay (Baltic Sea) in northern Germany. Similar variations of scintillation versus ASTD were observed and these measurements are presented in this paper.
In this work we reexamine the EOPACE data with the aim of proposing a bulk C2n model that matches the data more closely. Our primary goal, therefore, is to describe the optical refractive index structure parameter C2n using the MOS theory as a starting point, and then extending it by incorporating bulk parameters that are outside its scope. We emphasize that our focus is on bulk models that use average meteorological measurements to estimate quantities related to the turbulent fluxes (see the appendix). We do not consider models using high-frequency eddy covariance estimates of those fluxes. Through some assumptions on the causes of the discrepancies discovered by FDZB, changes are made to our current C2n expression, and a Bayesian method is applied to fit the new parameters in agreement with the observations. Section 2 reintroduces the MOS theory applied to the calculation of C2n in the marine surface layer and presents the basic elements of the Bayesian method. In section 3, a new expression is derived for C2n based on the EOPACE data and is tested against the observations made during VAMPIRA. We discuss the physical implications of our new model in section 4 and conclude in section 5.
a. Monin–Obukhov similarity theory and C2n
The perennial theory for describing turbulence in the atmospheric surface layer is MOS theory. The physical assumptions of MOS theory are that the surface layer is statistically horizontally homogeneous and stationary, which essentially implies that the fluxes of sensible heat, momentum, and humidity are constant throughout the surface layer. This assumption allows us to define characteristic scales of velocity u*, potential temperature θ*, and specific humidity q* (Kaimal and Finnigan 1994), given by
The primed quantities represent turbulent fluctuations of horizontal and vertical wind velocities (u and w, respectively), potential temperature (θ), specific humidity (q), and the overbars represent the average over the turbulent fluctuations.
We use the traditional Reynolds decomposition in which a random field R(x, t), defined over spatial coordinate x and time t, is split into average and fluctuating fields, R(x, t) = R(x, t) + R′(x, t) and = 0. The characteristic scales are constants over space and time in conventional MOS theory. However, as we wish to consider cases with variable fluxes, the scales in Eqs. (1)–(3) should be understood as variable scales that depend upon local fluxes. The structure function for the field R is defined by
where l is a displacement vector. If the turbulence is homogeneous, isotropic, and inertial, then the structure function depends only on the magnitude of l such that (Kaimal and Finnigan 1994)
where C2R is the structure parameter for the field R associated with its structure function.
A normalized height coordinate ξ = z/L (or stratification) is defined, where z is the height above ground and L is the Monin–Obukhov length. The normalized coordinate ξ is representative of the ratio of the buoyancy production of turbulent kinetic energy (TKE) to the shear production of TKE. General MOS theory states that the characteristic scale, along with the height and the normalized height, determines every profile associated with that scale. The structure parameters for potential temperature and specific humidity are given as
where the functions gθ, gq, and gθq are the similarity functions for the potential temperature, specific humidity, and cross structure parameters, respectively. Of these, it is the potential temperature structure function that has been studied the most (Garratt 1992). Its similarity function may be modeled according to Edson et al. (1991):
where C0 = 5.8, CU = 7, and CS = 2.4. Furthermore, it is commonly assumed that gθ = gq = gθq = g (see, e.g., Andreas 1988; Hill 1989). We follow this convention except that we fix gθq = rgθgq = rg, where r is the correlation coefficient between the potential temperature and specific humidity fluctuations. We set r = 0.8, as was done by some authors (Forand 1999; Frederickson et al. 2000) to better reflect measured values of that coefficient (Kohsiek 1982; Priestley and Hill 1985; Andreas 1987).
We need a formula connecting the structure parameters [Eqs. (6)–(8)] with the refractive index structure parameter C2n. For infrared wavelengths, the refractivity for moist air, Z = (n − 1) × 106, where n is the refractive index, can be given by
where T is the temperature in kelvins, Pa (hPa) is the partial pressure of dry air, Pw (hPa) is the partial pressure of water vapor, and A and B (K hPa−1) are wavelength-dependent coefficients [see Andreas (1988) and Beland (1993) for the detailed wavelength dependence]. We must now consider the fluctuations of Z. When doing this it is usually assumed that the total pressure does not fluctuate (i.e., there is instantaneous pressure equalization), such that P′ = P′a + P′w ≈ 0 (Andreas 1988) and the equation of state for water vapor yields
where Rw = 461.5 J kg−1 K−1 is the gas constant for water vapor, P is the total pressure, and ρw is the water vapor density. In terms of potential temperature and specific humidity, Eq. (11) becomes
where P0 = 1000 hPa is the reference pressure used in the definition of potential temperature, Rd = 287 J kg−1 K−1 is the gas constant of dry air, and κ = 0.286 and ɛ = 0.608 are dimensionless constants. The turbulent refractivity fluctuations can now be written to first order as
Since the refractivity structure function is proportional to the square of its turbulent fluctuation, it follows that the refractivity structure parameter has the same quadratic form, namely,
Equation (15) allows us to express the refractivity structure parameter in terms of MOS theory through the use of Eqs. (6)–(9). Since we will be using scintillation data from two different trials with different measurement heights, it is convenient to define a normalized structure parameter
Various implementations of the MOS theory have been presented in the literature for the computation of θ*, q*, and L. The various methods differ mainly in the selection of the stability functions and the values of coefficients and constants. Details of the calculations used in this paper are given in the appendix.
An alternative formulation of C2n developed by Tatarskii (1961) uses the mean refractive index profile
where LF is proportional to what Tatarskii calls the outer scale. Some authors (Walters and Kunkel 1981; Beland and Brown 1988) have taken LF to be proportional to the outer scale LO, where the structure function is roughly twice the refractive index variance C2nL2/3O ≈ 2σ2n. Beland and Brown (1988) questioned whether this is actually the case. In fact, Tatarskii may have had something else in mind. He arrives at Eq. (17) by considering random variations for the refractive index due to turbulence, δn2 = C2nl2/3, “systematic” variations due to changes in the average field, δn2 = (∇n)2l2 ≈ (∂n/∂z)2l2, and equating the two. The scale LF is therefore the scale where the random turbulent variation equals the systematic variation due to the average profile or, alternatively, the scale where it is impossible to distinguish turbulent fluctuations from changes in the average background field. Viewed this way, LF is the scale where we can no longer meaningfully measure turbulent fluctuations. Using the same reasoning as before, we can write
Then, MOS theory enables us to write the gradients of the mean profiles as functions of the corresponding characteristic scale and a scaling function
where k is the von Kármán constant. Assume that ϕθ = ϕq = ϕ, where
and combining this result with Eq. (16), we obtain what we call a Monin–Obukhov (MO) fluctuation length, LF = zΛ(ξ), which is a function of z and ξ and where
This function can be considered as the MO fluctuation length normalized with respect to the height (see Fig. 1). The Tatarskii formulation can provide an intuitive explanation for the discrepancies between the MO refractivity structure parameter and the values obtained from the transmission paths at EOPACE and VAMPIRA. Namely, the interaction between the atmospheric surface layer and the sea surface causes the MO fluctuation length to be shorter or longer than what it would be over land (the environment for which MOS theory was originally conceived). As will be seen later [Eq. (29)], the Tatarskii formulation justifies the use of a positive multiplicative function to characterize the scintillation deficit.
b. Bayesian methods
Bayesian methods are most useful when we try to estimate the value of an inaccessible parameter β from direct observations D, given that there exists a probabilistic relationship between them in the form of a conditional probability density (or likelihood) function ρ(D|β) (the probability of observing the value D given the value β). Bayesian methods also assume that we have information I about β that allows us to state our beliefs in the form of a probability density, ρ0(β|I), prior to considering any data (hence the term “a priori distribution”). After we obtain the data, Bayes’s rule allows us to update our a priori beliefs into a posteriori beliefs,
The a posteriori probability density, ρ1(β|D, I), therefore represents our belief about the parameter given our prior information and the newly observed data. It is this ability to introduce prior information that sets Bayesian methods apart from conventional methods (e.g., least squares and Chi square). There are many things we can do with the a posteriori probability density, such as finding the mean or standard deviation. In this work, we will simply find its mode, or the value of βmax that maximizes ρ1(β|D, I). If we have not one but a set of independent data D = (D1, . . . , Dn), then Eq. (23) generalizes to
If the number of observations n is large, then the likelihood function for the dataset (the product of the individual likelihood functions) is often very sharply peaked about its maximum likelihood value, βml. If the prior probability density is relatively flat around that value, then ρ0 can be treated effectively as a constant and removed from the integral in the denominator, thereby canceling itself in the numerator. Then, the maximum a posteriori probability value effectively equals the maximum likelihood value, βmax ≈ βml. The prior information about the parameter is thereby rendered irrelevant by the large number of observations, making the Bayesian approach superfluous. However, such a conclusion only holds when the prior probability is relatively smooth and continuous. If, for instance, there are good reasons to believe that the parameter must be positive, then this information can be incorporated by setting ρ0(β|I) = 0 for β < 0 such that no amount of data can force a negative result, something not all techniques can guarantee.
There is a certain amount of subjectivity in the choice of information about the parameter and the way it is included in the prior distribution (i.e., the specific shape of the prior) and into the results of a Bayesian analysis. However, Bayesians would argue that such subjectivity is always present but hidden, to give the analysis an aura of objectivity. It is better, they argue, to openly admit and quantify one’s prior beliefs about the estimated parameter and to allow others to judge the results (see D’Agostini 2003).
3. Data analysis
a. The EOPACE and VAMPIRA trials
In our analysis we consider scintillation data collected during two measurement campaigns: the U.S.-led EOPACE campaigns conducted in San Diego during the spring and autumn of 1998, and the North Atlantic Treaty Organization (NATO)-led VAMPIRA campaign conducted across Eckernförde Bay in northern Germany in the spring of 2004.
During EOPACE, the Space and Naval System Warfare Command Center, San Diego (SSC-San Diego) measured scintillation of a midwave infrared signal transmitted in a continuous manner across San Diego Bay. The transmitter was mounted 5 m above mean sea level (MSL) at the Naval Amphibious Base, Coronado, and the receiver was mounted at about the same height at the Naval Submarine Base, Point Loma, 7.0 km from the transmitter. The transmitter contained a 1200-K blackbody at the focal plane of an F/6 Newtonian telescope, projecting a conical beam with a full angle of 5.4 mrad. Details on the transmitter–receiver setup are given by FDZB.
Meteorological measurements required for bulk model computations were obtained simultaneously using an instrumented buoy deployed by the Naval Postgraduate School (NPS) in Monterey at about the midpath between the transmitter and the receiver. Wind speed and direction were measured on the buoy by a sonic anemometer 4.9 m MSL and the air temperature and humidity were measured 3.6 m MSL. The bulk sea temperature was measured with a thermistor placed within a float attached alongside the buoy’s hull. For this analysis, 10-min averages of the 1-Hz raw data are used for the bulk computations of C2n.
During VAMPIRA, Defense Research & Development Canada (DRDC)—Valcartier investigated scintillation by obtaining high-speed visible imagery of light sources using a high-speed digital camera. The camera was mounted on a pier in Surendorf, Germany, 6.6 m MSL. Halogen light sources were installed 8.2 km across Eckernförde Bay in Booknis Eck, 8.5 m MSL. Air temperature and humidity along with sea temperature were measured concurrently in the middle of the Bay using a buoy from the Netherlands Organization for Applied Scientific Research (TNO-FEL). The air temperature and humidity probe was 5.2 m MSL and the sea temperature was measured about 1 m below the surface. Wind speed and direction were measured about 5.5 m MSL near the camera system [see Heemskerk (2005) and Seiffer and Stein (2005) for more information on VAMPIRA].
The high-speed camera could acquire images at rates ( fC) between 50 and 2000 frames per second and 2000 image sequences were recorded approximately every half-hour. As turbulent phenomena such as scintillation, image displacement, and blurring occur with characteristic frequencies that are proportional to the wind component perpendicular to the propagation path, wind measurements from an anemometer near the camera were used to determine the transverse wind speed and to set the camera’s frame rate. The frame rate was set to at least 10 times the scintillation characteristic frequency ( f0), and the acquisition period (2000/fC) was set to at least 10 times the characteristic period (1/f0). This ensured that each sequence adequately captured the full spectrum of turbulent imagery phenomena (no aliasing while resolving the low frequencies). For the postprocessing, a special DRDC program was used to remove the background surrounding each light, and then calculate the integrated intensity and center-of-mass for each light, for each image of every sequence (see Potvin et al. 2007). From these, estimates of the average and variance of the integrated intensity of the lights for each sequence provided the scintillation measurements from which C2n values were inferred. The center-of-mass variances, which are also proportional to the C2n, were also used to corroborate the C2n measurements.
b. EOPACE scintillation data
Figure 2 presents the normalized structure parameter S derived from the EOPACE scintillation data along with bulk estimates obtained using the model described in section 2. The model assumes identical similarity functions, gθ = gq = gθq, using coefficients proposed by Edson et al. (1991) and is similar to the bulk model used by FDZB. Two problems are immediately apparent, as noted by FDZB. First, near neutrality the model predictions drop far below the observations. Second, the model clearly overestimates the observations when the surface layer is strongly unstable and, especially, when it is appreciably stable. FDZB suggests various possible explanations for the observed discrepancies. Assumptions made in this paper for the derivation of a new expression for C2n are compatible with FDZB explanations.
In this paper, two changes are made to the current C2n model to match the observations over water more closely. First, the fluctuation length will be allowed to vary and, second, assuming that the discrepancies near neutrality are mainly due to nonuniform fluxes along the propagation path, variations of (or uncertainties about) θ* and q* will also be taken into account.
c. Flux nonuniformity
An obvious physical effect perturbing the modeling of C2n is the fact that the meteorological conditions are not uniform along the propagation path. This is particularly important near neutrality where the sensible heat flux is expected to drop significantly along with C2n. Expectedly, even small variations in the surface fluxes along the path can prevent the C2n from falling along with the bulk model. If the refractive index structure parameter varies along the propagation path, C2n (x), then the scintillometer will measure a path-weighted structure parameter
where W = [x(LP − x)]5/6 is the scintillometer’s pathweighting function and LP is the pathlength. Similarly, the characteristic scales can also be functions of the propagation path, θ*(x) and q*(x), where x goes from 0 to LP, their path-averaged values θ̃* and q̃* are determined by
and their path variances are given by σ2θ = − θ̃2* and σ2q = − q̃2*. For simplicity, we shall assume that = 0 such that the normalized path-weighted structure parameter becomes
which can be decomposed as S̃ = S0 + SV, where
In strict terms, the variation in the fluxes should also cause variations in the normalized height ξ. However, we neglect this, as the similarity functions are not very sensitive to such variations. Furthermore, from here on we will drop the tildes so that S, θ*, and q* should be understood as S̃, θ̃*, and q̃*, respectively. It is important to note that at neutral stability there may be temperature and humidity fluctuations that are not related to flux variations, such as advection of the temperature variance along the path. It is possible that the path variances estimated from the data include such effects.
d. Anomalous deficit for large ASTDs
The expression of C2n needs a further correction to describe the observed deficit for large ASTDs. We can exclude the attenuation of the scintillation due to multiple scattering of the radiation from the turbulent refractive index fluctuations, as FDZB used a large aperture scintillometer, resistant to such effects. According to the theory of Wang et al. (1978) this only begins at a normalized structure parameter value of 0.144 and the EOPACE scintillation data of Fig. 2 is well below this threshold. As will be shown in the following section, a similar deficit was observed at VAMPIRA (see Fig. 3), and the deduced C2n proved to be in agreement with the observed image displacements. The physical cause for the deficit is not easy to determine. In this paper we will attempt to account for it in a purely empirical manner, within the context of a variable fluctuation length hypothesis.
The anomalous deficit depends most strongly on θυ*, the virtual potential temperature characteristic scale. No obvious correlation has been found between the deficit and the stratification, ξ, which suggests that this phenomenon lies outside MOS theory’s domain. The fact that only points for which |ξ| ≤ 1 are used in the analysis, and therefore for which MOS theory should apply, reinforces this conclusion. We model changes in the fluctuation length by a correction function F(θυ*) such that L′F = exp[F(θυ*)]LF(ξ). Recalling that the structure parameter is proportional to the fluctuation length to the 4/3 power [see Eq. (17)], our empirical model will have the form
e. Prior probability and Bayesian analysis
The problem can now be formulated in Bayesian terms, as introduced in section 2b. The first task is to define a functional form for the function F. Because of the considerable scatter of the data, we can plausibly use a wide variety of functions to describe F. Furthermore, as the points are not evenly distributed along the abscissa, this means that the points concentrated about the origin will dominate the fit at the expense of the points at the extremities. Worse, the quality of the mismatch at the extremities is very sensitive to the mathematical form chosen for the function F. To avoid these difficulties, we begin by modeling F (see Fig. 6) as a series of constants spanning certain intervals in θυ*, where F = f1 over the interval − 0.2 < θυ* < −0.1, F = f2 for −0.1 < θυ* < −0.05, F = f3 for −0.05 < θυ* < 0, and similarly for θυ* positive to give us six constants (see Fig. 6). By breaking up F into independent segments, we assure a uniform quality of fit over the entire span of θυ*, regardless of the nonuniform distribution of points.
Once the constants have been estimated, we then find representative θυ* values for each interval by finding the average of these values within the interval. Subsequently adding the two extreme points, (−0.2, f1) and (0.2, f6), gives us two 8D vectors from which we generate a continuous correction function using a natural cubic spline interpolation with the first-order derivatives set to zero at the end points. Since we know nothing about the anomalous deficit, we postulate a uniform prior distribution for its coefficients, ρ0( f1, . . . , f6|∅︀) = 1, over all values of the coefficients. Strictly speaking, such a prior is improper as its integration yields infinity (and not unity). Nevertheless, improper priors are often easier to treat mathematically and nonetheless, as in this case, lead to proper posterior distributions. This can be understood by looking at Eq. (23). If the prior is constant, then it can be removed from the integral in the denominator and is cancelled by the numerator. If the integral of the likelihood function with respect to the parameter in the denominator is finite, then the resulting posterior distribution integrates to unity and is proper.
These deficit coefficients are assumed to be “universal” constants in that a single value applies for the entire dataset. The path variances are a different matter, as they refer to the prevailing meteorological conditions for a particular datum; thus, we expect to have a different path variance value for each scintillation measurement. We must then find a prior distribution for path deviation vectors σθ = (σθ,1, . . . , σθ,n) and σq = (σq,1, . . . , σq,n), where n (very large for EOPACE) is the number of points in the dataset.
However, rather than finding an optimal solution in an n + 6 dimensional space, we will instead reduce the path deviation vectors to their sufficient statistics. In other words, let us suppose that the path deviation vectors have a prior distribution of the form
which can be rewritten as
where the angle brackets represent the average over the entire dataset (〈σ2θ〉 = n−1Σni=1σ2θ,i and so on). The statistics 〈σ2θ〉, 〈σθσq〉, and 〈σ2q〉 are sufficient in the sense that they determine the probability for the entire set of path deviations.
We can further simplify Eq. (31) by the decompositions 〈σ2θ〉 = 〈σθ〉2 + χ2θ, 〈σθσq〉 = 〈σθ〉〈σq〉 + χθq, and 〈σ2q〉 = 〈σq〉2 + χ2q. From Eq. (31) it is clear that we can also decompose the prior such that ρ0(σθ, σq) = ρn0(〈σθ〉, 〈σq〉)ρn0(χ2θ, χθq, χ2q). We can achieve further simplification by assuming that the second-order statistics χ2θ, χθq, and χ2q are very small compared to the averages 〈σθ〉 and 〈σq〉, so that only the averages play a role in modeling C2n, and enabling us to integrate out the second-order statistics. Thus, we have reduced a problem involving 2n parameters to one with two parameters (〈σθ〉 and 〈σq〉) with a prior distribution ρ0(〈σθ〉, 〈σq〉) such that
The variables μθ, μq, and η are parameters that determine the shape of the prior, which is itself a distribution of parameters. We therefore call these variables hyperparameters (Gelman et al. 2004), which we must determine based on some a priori knowledge of the nonuniformity. During previous measurements, SSC-San Diego had made nine runs along the propagation path in San Diego Bay using an instrumented Boston Whaler, measuring the ASTD and air–sea specific humidity difference (ASQD) along the way (Fig. 4 shows a sample). The overall standard deviation for the ASTD was about 0.59°C and for the ASQD about 4.12 × 10−4 kg kg−1. This translates roughly as standard deviations of 2.05 × 10−2 K for the potential temperature scale, and of 1.65 × 10−5 kg kg−1 for the specific humidity scale. Given the errors associated with the boat measurements, we take these values as maximums for the possible values of 〈σθ〉 and 〈σq〉. However, even without measurement errors these estimates would be too high, as they do not take into account the “footprint” of the fluxes. The flux at a certain location and a certain height above the surface is the integration of the surface flux over a certain surface area or footprint (Schmid 1997). This means that the ASTD and ASQD values shown in Fig. 4 should be smoothed to reflect the footprint effect, thereby reducing the estimated standard deviations. Consequently, we would be very skeptical of any result that showed that the average path deviation was greater than the boat value, and this attitude should be reflected in our prior.
We also know that the Bowen ratio (ratio of sensible heat flux over latent heat flux) is about 0.2 over water (Garratt 1992) and is confirmed by the EOPACE data, which give a median value of about 0.2. Since the sensible heat flux is proportional to θ* and the latent heat flux to q*, the Bowen ratio imposes proportionality between these scales and between their path deviations. To respect the Bowen ratio and not exceed the limits set by the boat measurements, we set μθ = 8.55 × 10−3 K and μq = 1.65 × 10−5 kg kg−1. Furthermore, the EOPACE dataset reveals a correlation of about 0.85 between θ* and q*, and so from this we simply set η = 0.85. The estimates of μθ, μq, and η are admittedly rough. However, it is worth noting that the end result does not depend strongly on these hyperparameters, so their exact value is not crucial.
Now that we have found an acceptable form for the prior distribution, we must find the likelihood function ρ(SD|S), where SD is the normalized structure parameter derived from the EOPACE scintillation data. It must satisfy several requirements; namely, it should reflect that SD can only be positive, and that it has a wide dynamic range [so the likelihood function for ΩD = log(SD) should have an acceptable form], and it should be unbiased, such that ∫SDρ(SD|S)dSD = S. Since we will be dealing with the logarithm of SD, and will find the parameters with the highest a posteriori probability, we would like the likelihood for ΩD
where Ω = log(S), to have a maximum-likelihood value Ω̃D such that Ω̃D = Ω.
A good choice that satisfies all these requirements is the gamma distribution,
where m is a parameter determining the variance of the distribution, and the logarithmic version has the form
Figure 5 shows that m = 9 is a good choice, yielding a function that is reasonably Gaussian-looking. It also shows a likelihood function estimated from the EOPACE scintillation data and the modified bulk model using the correction function shown in Fig. 6. There is reasonable agreement between the two functions. Although the estimated likelihood function is more peaked than the chosen gamma distribution, this does not significantly affect the Bayesian analysis. With our prior distribution and the likelihood function of Eq. (35), we can use Eq. (24) to find the eight-dimensional posterior probability density function ρ(〈σθ〉, 〈σq〉, f1, . . . , f6|ΩD, I).
Using a numerical maximization procedure, we find the parameter values that maximize this function (see Table 1). Using these parameter values in our model, Eq. (29), gives our optimal model estimates for the normalized structure parameter. Figure 7 shows the new model calculations along with the EOPACE values. As expected, there is a marked improvement with respect to the conventional model calculations shown in Fig. 1. The improvement is also visible in the double scatterplot of Fig. 8, where the modified bulk model is much better aligned with respect to the diagonal than the conventional bulk model. The improvement is also evident as a function of the virtual potential temperature scale, shown in Fig. 9.
With the conventional bulk model, we see that the median scintillation deficit is steady at about −6 dB for strong positive values of θυ*, whereas for negative values the scintillation deficit seems to get progressively worse, with a minimum at about −4.5 dB. The peak close to θυ* = 0 is the conventional bulk model dipping well below measured scintillation values at neutrality. The modified model on the whole stays close to the 0-dB line. It falls to a minimum of about −2.5 dB at θυ* = −0.05 and at θυ* = 0.03, but in both instances the 90th percentile intervals cover the 0-dB line. The only exception to this is at θυ* = −0.175, where the scintillation deficit falls to about −2 dB.
f. Application of the modified bulk model to VAMPIRA
Figure 10 shows the modified bulk model applied to the VAMPIRA data, where we show all 100 scintillation data points and all of the meteorologically derived data (874 points) to compensate for the relative scarcity of data. We use the same anomalous deficit function as for the EOPACE data, but have performed a Bayesian regression and obtained slightly different optimal path deviations. The VAMPIRA data in Figs. 3 and 10 are displayed using the same scale as the EOPACE data in Figs. 2 and 7 for easy comparison. We notice an improvement with respect to the conventional bulk model, shown in Fig. 3, for stable conditions. Such improvement does not appear under unstable conditions except maybe at larger temperature scales. However, the small number of points prevents us from coming to any firm conclusions.
It is worth mentioning that the ASTD had a smaller range at VAMPIRA (from −2.6° to 2.6°C) than at EOPACE (from −3.8° to 8.9°C), and VAMPIRA was slightly windier on the whole (wind speeds ranging from 0.5 to 12.4 m s−1) than EOPACE (0.3–9.0 m s−1). Furthermore, the winds for the VAMPIRA data tend to be well clustered around moderate speeds and to come from the open sea for unstable conditions and to be evenly distributed over a wide range of speeds and directions for stable conditions. This suggests that the sampling of meteorological conditions is very different between unstable and stable conditions. By contrast, the EOPACE dataset has fairly similar wind sampling for both the unstable and stable conditions (generally from the land with low to moderate wind speeds). To these meteorological considerations, we must add subtle geographical considerations (the form of the respective bays and their effects on wind flow, thermal properties, etc.) when comparing these datasets.
The first thing to point out is that the virtual potential temperature scale in Fig. 3 is much less variable for the VAMPIRA scintillation data than for the EOPACE data. This is largely due to the difference in height between the setups, as only points with an absolute normalized height coordinate less than unity are shown, which leads to a more restricted range of virtual potential temperature scales.
The height difference between the EOPACE and VAMPIRA measurements might also explain why the VAMPIRA data seems to have a more pronounced V shape than the EOPACE data. Another plausible explanation is that under neutral conditions at VAMPIRA the winds are systematically stronger (about 8 versus 3 m s−1) and blow from the sea, which could lead to greater homogeneity, translating into reduced path variances in our model. Figure 11 shows the same plot as in Fig. 10, but using significantly lower path deviations and for which better agreement is seen.
From Figs. 3, 10 and 11, one may be tempted to conclude that, for moderately unstable conditions, both the conventional and modified bulk model perform much less well at VAMPIRA than at EOPACE. This may suggest that the modeling of the anomalous deficit under stable conditions for the EOPACE data can carry over to the VAMPIRA data, but does not under unstable conditions. This can be alleviated somewhat by adjusting the path deviations (see Fig. 11). It is also worth noting that the EOPACE data (Figs. 2, 7) exhibit variability over about one decade and that for a given temperature scale at EOPACE, under unstable conditions, the scintillation indices were found to be correlated with the wind speed; in general the greater the wind speed, the greater the scintillation. At VAMPIRA, all measurements under unstable conditions correspond to an average wind speed of 4 m s−1. It is interesting to note that for this wind speed, the EOPACE measurements are less than the model predictions by about one-third of a decade, just like at VAMPIRA.
a. Physical origin of the anomalous deficit
As mentioned previously, the anomalous deficit is not an attenuation of the scintillation caused by multiple scattering, and as we limited ourselves to points |ξ| < 1, it is not because we applied MOS theory beyond its domain of validity. This last statement should be qualified because MOS theory was applied in a coastal region, with all its inherent nonuniformities that are antithetical to the very premises of MOS theory. Nevertheless, we believe the main factor is that the observations were done over water. The sea surface is not static and maintains a dynamic relationship with the air above it, unlike the mostly immovable surface of dry land. The scintillation index, from which we obtain estimates of C2n, depends on an integral over the entire turbulent refractive index spectrum (Wang et al. 1978), the form of which should depend in some way on the fluctuation length.
We can make the preceding discussion more quantitative by following the analysis of Andreas (1987) for the potential temperature structure parameter C2θ. This parameter depends on the energy dissipation ɛ, related to the TKE, and the dissipation rate of the potential temperature variance Nθ as C2θ ∝ Nθɛ−1/3. The MOS theory assumes that there is no advection of TKE and no wave interaction [see Chang and Cheng (1972) for an analysis of the turbulent kinetic energy budget that includes wave interaction terms]. Assuming a stationary surface layer and that wave-induced perturbations are independent of turbulent fluctuations, the energy dissipation can be written as
where ɛ0 is the MOS energy dissipation, Adv is the advection term, and WI is the wave interaction term. Performing the same treatment to the temperature variance dissipation, we may write
where Nθ0 is the MOS potential temperature variance dissipation. If we suppose that the conventional bulk model gives only C2θ0 ∝ Nθ0ɛ−1/30 and the scintillation measures C2θ, then we have an approximate account of the physical origin of the anomalous deficit. But is this really what the conventional bulk model gives us? It deduces everything from the average profiles of wind and temperature. Could the average profiles be affected by advection and wave interaction? Such questions are beyond the scope of this study. It is worth pointing out that the advection and wave interaction terms could be active under neutral conditions and could compete with the flux nonuniformity in explaining the absence of the large dip under those conditions.
We still wonder about the functional dependencies of the correction function. For instance, given the wave interaction terms in Eqs. (36) and (37) it seems reasonable to assume a dependence on the sea state, but which parameters? If the sea state is involved, then it may depend on a characteristic wave height, wavelength, or period. The deficit may depend upon quantities related to the nonuniformity of the environment, such as the advection terms of TKE and potential temperature variance in Eqs. (36) and (37) or specific humidity variance advection from the land onto the sea. Whatever the case may be, there is no reason to expect the advection and wave interaction terms to conform to Monin–Obukhov similarity. Clearly more study is needed involving a larger number of environmental quantities than hitherto thought necessary.
b. Generalizing the results
A modified bulk model was empirically derived using the abundant EOPACE data, and then applied to the VAMPIRA data to see how well the modified model performed. As discussed above, it is likely that in general the path deviation parameters are not transportable from one dataset to another, as it would logically depend on the setup (e.g., height of observations) and local conditions (geographical and meteorological), thereby making this part of the model site-specific. The physical cause of the deficit is yet to be clarified. Nevertheless, it is interesting to note that a deficit is observed for both EOPACE and VAMPIRA and that an appreciable improvement in the bulk predictions at VAMPIRA can be obtained for stable conditions by applying the correction function F derived from EOPACE data.
In this work, we created a modified bulk model for C2n using the EOPACE dataset and Bayesian methods, and then applied this model to predicting the C2n values from the VAMPIRA dataset. The modified model includes the effects of nonuniform fluxes and accounts for an anomalous scintillation deficit for large ASTDs. It can account for both datasets better than conventional bulk models. The anomalous deficit suggests how a model for C2n must extend beyond MOS theory for the above-water coastal environment. In our model, the deficit has been made dependent upon the virtual potential temperature scale, the only apparent dependent factor in our analysis. Although we postulate a connection with a characteristic scale of the turbulent spectrum, we nevertheless do not know the exact physical nature of this deficit; consequently, we have yet to determine the entire set of meteorological parameters controlling its behavior. This should be the focus of future work in this domain.
The authors gratefully acknowledge the efforts of all those who made the EOPACE and VAMPIRA datasets possible. In particular, regarding EOPACE, the authors thank Dr. Carl Zeisse from the SSC-San Diego for the EOPACE scintillation measurements, Dr. Stephen Hammel for his assistance in acquiring the data, and Dr. Douglas Jensen for his organization of all the campaigns. With regard to VAMPIRA, the authors thank Dr. Karin Stein and the TG32 NATO group for the organization and coordination of the campaign. Thanks are also given to Luc Gauthier of DRDC—Valcartier, for his tireless efforts during Valcartier’s efforts during both EOPACE and VAMPIRA. The data analysis presented in this paper was funded by Defence Research & Development Canada of the Department of National Defence and by the Direction Générale de l’Armement in France.
Calculation of the Scaling Parameters and the Monin–Obukhov Length
The estimation of C2n using Eq. (15) requires prior knowledge of the scaling parameters θ*, q*, and L— the Monin–Obukhov length (also known as stability length). The inputs (hereinafter called bulk parameters) needed for the calculation are the atmospheric pressure P1, the air temperature T1, the relative humidity H1 measured at a height z1 above the sea surface, the sea surface temperature T0, and the wind speed U2, measured at a height z2. (For notation in the appendix, subscript 0 refers to the parameters taken at the sea surface, whereas subscript n refers to a value at height zn.) Moreover, it is also assumed that the wind speed at the sea surface is zero (U0 = 0), and the relative humidity at the sea surface is given by
where s is the salinity expressed in grams per kilogram.
The relative humidity values also need to be converted into specific humidity. To this end, we recall that the relative humidity is defined as
where Pw and Pws are the water vapor pressure and the saturated water vapor pressure, respectively. As discussed by Forand (1999, 2007), the saturated vapor pressure (hPa) can be expressed using a simple function of the temperature T (K):
The specific humidity q (kg kg−1) depends on the total pressure P (hPa), and the water vapor pressure as given by
Finally, as within the marine surface boundary layer, the vertical profile of the total pressure is quite linear and can be expressed as
where g (≈9.807 m s−2) is the gravitational acceleration and Rd (≈287.05 J kg−1 K−1) is the gas constant for dry air.
The Monin–Obukhov length L is expressed as
where u* is the friction velocity, k is the von Kármán constant, θυ* = θ* + 0.61θq* is the virtual potential temperature scaling parameter, and θ and θυ are the average potential temperature and average virtual potential temperature, respectively, which are estimated by θ = (θ1 + θ2)/2 and θυ = (θυ,0 + θυ,1)/2. As discussed in Frenzen and Vogel (1995), many studies tend to show that k actually varies depending on the surface characteristics (more specifically on the thickness of the laminar sublayer). In our model k is taken to be 0.41.
The MOS theory stipulates that the mean vertical profiles of wind speed (U), potential temperature (θ), and specific humidity (q) are given by
where z0u, z0t, and z0q (called roughness lengths) are the heights at which the mean vertical profiles of U, θ, and q reach their surface values. The profile stability functions Ψ have different empirical expressions for the unstable case (L < 0) or stable case (L > 0). For unstable conditions, we use the functions derived by Edson et al. (1991):
whereas for the stable cases, the function proposed by Kondo (1975) is used:
Unlike the other concurrent functions proposed in the literature for stable conditions, the Kondo functions are attractive since they are logarithmic and thus limit the increase of the profile with height. These stability functions were derived empirically from measurements over land (with k equal to 0.4). The validity of these expressions to describe optical effects over the sea has been of concern for scientists and engineers [an assessment of their performance for predicting refraction effects is presented in Dion et al. (2005)]. Results obtained so far suggest that adjusting k to a value of 0.41 leads to improved model predictions.
Through a parameterization of the roughness lengths (see below), and by evaluating the three basic Eqs. (A7)–(A9) at the measurement heights, we obtain in conjunction with Eq. (A6) a set of four equations for the four unknown parameters: u*, θ*, q*, and L. The solution is obtained through an iterative process.
The momentum roughness length z0u can be expressed according to Fairall et al. (1996) as
where αc is Charnock’s “constant” and ν is the kinematic viscosity of air. As suggested by Smith et al. (1992), we assume that αc = 0.011 for open-sea conditions and αc ≈ 0.017 in coastal regions (the value used with the EOPACE and VAMPIRA data). Numerical values listed by Garratt (1992) show that the kinematic viscosity can be empirically expressed as a linear function of tc (expressed in degrees Celsius; thus tc = T − 273.15):
A simple way to determine the roughness lengths z0t and z0q is by making use of the neutral transfer coefficients CHN and CEN:
As mentioned by Garratt (1992), the experimental observations made under various wind conditions lead to the following values: CHN = CEN ≈ 1.1 × 10−3 (±15%) implying in particular that z0t = z0q. The dependence of z0u with the friction velocity u* as given by Eq. (A13) indicates that the calculation of the roughness lengths has to be included within the iterative process used to evaluate the scaling parameters.
Corresponding author address: Guy Potvin, Defence R & D Canada–Valcartier, 2459 Pie-XI Blvd., N. Québec, QC G3J 1X5, Canada. Email: firstname.lastname@example.org