## Abstract

This article examines some of the difficulties associated with the determination of *C*^{2}_{n} 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 *C*^{2}_{n} are discussed.

## 1. Introduction

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 *C*^{2}_{n}, 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 *C*^{2}_{n} obtained from standard environmental measurements.

Bulk models of *C*^{2}_{n} 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 *C*^{2}_{n} model that matches the data more closely. Our primary goal, therefore, is to describe the optical refractive index structure parameter *C*^{2}_{n} 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 *C*^{2}_{n} 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 *C*^{2}_{n} in the marine surface layer and presents the basic elements of the Bayesian method. In section 3, a new expression is derived for *C*^{2}_{n} 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.

## 2. Theory

### a. Monin–Obukhov similarity theory and C^{2}_{n}

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 *C*^{2}_{R} 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*_{θ}, *g*_{q}, 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 *C*_{0} = 5.8, *C*_{U} = 7, and *C*_{S} = 2.4. Furthermore, it is commonly assumed that *g*_{θ} = *g*_{q} = *g*_{θq} = *g* (see, e.g., Andreas 1988; Hill 1989). We follow this convention except that we fix *g*_{θq} = *r**g _{θ}g_{q}* =

*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 *C*^{2}_{n}. For infrared wavelengths, the refractivity for moist air, *Z* = (*n* − 1) × 10^{6}, where *n* is the refractive index, can be given by

where *T* is the temperature in kelvins, *P _{a}* (hPa) is the partial pressure of dry air,

*P*(hPa) is the partial pressure of water vapor, and

_{w}*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 *R _{w}* = 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 *P*_{0} = 1000 hPa is the reference pressure used in the definition of potential temperature, *R _{d}* = 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

where

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 *C*^{2}_{n} developed by Tatarskii (1961) uses the mean refractive index profile

where *L _{F}* is proportional to what Tatarskii calls the outer scale. Some authors (Walters and Kunkel 1981; Beland and Brown 1988) have taken

*L*to be proportional to the outer scale

_{F}*L*, where the structure function is roughly twice the refractive index variance

_{O}*C*

^{2}

_{n}

*L*

^{2/3}

_{O}≈ 2

*σ*

^{2}

_{n}. 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,

*δ*

*n*

^{2}=

*C*

^{2}

_{n}

*l*

^{2/3}, “systematic” variations due to changes in the average field,

*δ*

*n*

^{2}= (

**∇**

*n*)

^{2}

*l*

^{2}≈ (∂

*n*/∂

*z*)

^{2}

*l*

^{2}, and equating the two. The scale

*L*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,

_{F}*L*is the scale where we can no longer meaningfully measure turbulent fluctuations. Using the same reasoning as before, we can write

_{F}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 where *D _{U}* = 20 (Edson et al. 1991) and

*D*= 6 (Kondo 1975). Placing this in Eq. (17) we obtain

_{S}and combining this result with Eq. (16), we obtain what we call a Monin–Obukhov (MO) fluctuation length, *L*_{F} = *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** = (*D*_{1}, . . . , *D*_{n}), 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 *C*^{2}_{n}.

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 ( *f _{C}*) 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 (

*f*

_{0}), and the acquisition period (2000/

*f*) was set to at least 10 times the characteristic period (1/

_{C}*f*

_{0}). 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

*C*

^{2}

_{n}values were inferred. The center-of-mass variances, which are also proportional to the

*C*

^{2}

_{n}, were also used to corroborate the

*C*

^{2}

_{n}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*_{θ} = *g*_{q} = *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 *C*^{2}_{n} are compatible with FDZB explanations.

In this paper, two changes are made to the current *C*^{2}_{n} 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 *C*^{2}_{n} 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 *C*^{2}_{n}. Expectedly, even small variations in the surface fluxes along the path can prevent the *C*^{2}_{n} from falling along with the bulk model. If the refractive index structure parameter varies along the propagation path, *C*^{2}_{n} (*x*), then the scintillometer will measure a path-weighted structure parameter

where *W* = [*x*(*L*_{P} − *x*)]^{5/6} is the scintillometer’s pathweighting function and *L _{P}* 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

*L*, their path-averaged values

_{P}*θ̃*

_{*}and

*q̃*

_{*}are determined by

and their path variances are given by *σ*^{2}_{θ} = − *θ̃*^{2}_{*} and *σ*^{2}_{q} = − *q̃*^{2}_{*}. For simplicity, we shall assume that = 0 such that the normalized path-weighted structure parameter becomes

which can be decomposed as *S̃* = *S*_{0} + *S*_{V}, 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 *C*^{2}_{n} 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 *C*^{2}_{n} 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*(*θ*_{υ}_{*})]*L*_{F}(*ξ*). 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* = *f*_{1} over the interval − 0.2 < *θ*_{υ}_{*} < −0.1, *F* = *f*_{2} for −0.1 < *θ*_{υ}_{*} < −0.05, *F* = *f*_{3} 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, *f*_{1}) and (0.2, *f*_{6}), 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}( *f*_{1}, . . . , *f*_{6}|∅︀) = 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}Σ^{n}_{i=1}*σ*^{2}_{θ,i} and so on). The statistics 〈*σ*^{2}_{θ}〉, 〈*σ*_{θ}*σ*_{q}〉, and 〈*σ*^{2}_{q}〉 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 〈*σ*^{2}_{q}〉 = 〈*σ*_{q}〉^{2} + *χ*^{2}_{q}. From Eq. (31) it is clear that we can also decompose the prior such that *ρ*_{0}(*σ*_{θ}, *σ*_{q}) = *ρ*^{n}_{0}(〈*σ*_{θ}〉, 〈*σ*_{q}〉)*ρ*^{n}_{0}(*χ*^{2}_{θ}, *χ _{θq}*,

*χ*

^{2}

_{q}). We can achieve further simplification by assuming that the second-order statistics

*χ*

^{2}

_{θ},

*χ*, and

_{θq}*χ*

^{2}

_{q}are very small compared to the averages 〈

*σ*

_{θ}〉 and 〈

*σ*

_{q}〉, so that only the averages play a role in modeling

*C*

^{2}

_{n}, and enabling us to integrate out the second-order statistics. Thus, we have reduced a problem involving 2

*n*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 *ρ*(*S*_{D}|*S*), where *S _{D}* is the normalized structure parameter derived from the EOPACE scintillation data. It must satisfy several requirements; namely, it should reflect that

*S*can only be positive, and that it has a wide dynamic range [so the likelihood function for Ω

_{D}_{D}= log(

*S*

_{D}) should have an acceptable form], and it should be unbiased, such that ∫

*S*

_{D}

*ρ*(

*S*

_{D}|

*S*)

*d*

*S*

_{D}=

*S*. Since we will be dealing with the logarithm of

*S*, and will find the parameters with the highest a posteriori probability, we would like the likelihood for Ω

_{D}_{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}〉, *f*_{1}, . . . , *f*_{6}|**Ω**_{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.

## 4. Discussion

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 *C*^{2}* _{n}*, 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 *C*^{2}_{θ}. This parameter depends on the energy dissipation ɛ, related to the TKE, and the dissipation rate of the potential temperature variance *N*_{θ} as *C*^{2}_{θ} ∝ *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 *C*^{2}_{θ0} ∝ *N*_{θ0}ɛ^{−1/3}_{0} and the scintillation measures *C*^{2}_{θ}, 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.

## 5. Conclusions

In this work, we created a modified bulk model for *C*^{2}_{n} using the EOPACE dataset and Bayesian methods, and then applied this model to predicting the *C*^{2}_{n} 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 *C*^{2}_{n} 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.

## Acknowledgments

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.

## REFERENCES

**,**

^{2}

_{n}over snow and sea ice from meteorological data.

**,**

**,**

**,**

^{2}

_{n}en atmosphère marine. NASA Tech. Rep. AGARD-CP-567, 29-1–29-11. [Available from NASA Center for Aerospace Information, 800 Elkridge Landing Rd., Linthicum Heights, MD 21090-2934.]

**,**

**,**

**,**

**,**

_{n}

^{2}) over the ocean using bulk methods.

**,**

**,**

**,**

**,**

**,**

^{2}

_{T}, C

^{2}

_{Q}, and C

_{TQ}in the unstable surface layer, and relations to the vertical fluxes of heat moisture.

**,**

**,**

**,**

**,**

**,**

_{0}.

**,**

^{2}

_{n}.

**,**

### APPENDIX

#### Calculation of the Scaling Parameters and the Monin–Obukhov Length

The estimation of *C*^{2}_{n} 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 *P*_{1}, the air temperature *T*_{1}, the relative humidity *H*_{1} measured at a height *z*_{1} above the sea surface, the sea surface temperature *T*_{0}, and the wind speed *U*_{2}, measured at a height *z*_{2}. (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 *z*_{n}.) Moreover, it is also assumed that the wind speed at the sea surface is zero (*U*_{0} = 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 *P _{w}* and

*P*

_{ws}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 *R _{d}* (≈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

and

where *z*_{0u}, *z*_{0t}, and *z*_{0q} (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 *z*_{0u} 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 *t*_{c} (expressed in degrees Celsius; thus *t*_{c} = *T* − 273.15):

A simple way to determine the roughness lengths *z*_{0t} and *z*_{0q} is by making use of the neutral transfer coefficients *C*_{HN} and *C*_{EN}:

As mentioned by Garratt (1992), the experimental observations made under various wind conditions lead to the following values: *C*_{HN} = *C*_{EN} ≈ 1.1 × 10^{−3} (±15%) implying in particular that *z*_{0t} = *z*_{0q}. The dependence of *z*_{0u} 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.

## Footnotes

*Corresponding author address:* Guy Potvin, Defence R & D Canada–Valcartier, 2459 Pie-XI Blvd., N. Québec, QC G3J 1X5, Canada. Email: guy.potvin@drdc-rddc.gc.ca