## Abstract

Atmospheric turbulence is a primary concern for astronomers. Turbulence causes amplitude and phase fluctuations in electromagnetic waves propagating through the atmosphere, constraining the maximum telescope resolution and resulting in telescope image degradation. Astronomical parameters that quantify these effects are generically referred to as *seeing*. Adaptive optics (AO) is used to reduce image degradation associated with optical turbulence. However, to optimize AO, knowledge of the vertical profile of turbulence and overall (integrated) seeing is needed. In this paper, an optical turbulence algorithm is described that makes use of the information on turbulence kinetic energy provided by a planetary boundary layer scheme available in the fifth-generation Pennsylvania State University–NCAR Mesoscale Model (MM5). Optical turbulence data collected on Mauna Kea during the 2002 site monitoring campaign are used to validate the algorithm, which has been implemented in operational runs of MM5 at the Mauna Kea Weather Center.

## 1. Introduction

Atmospheric turbulence is a key challenge in ground-based astronomy because it dramatically impacts the angular resolution of a telescope. Small-scale temperature and moisture fluctuations in the atmosphere result in fluctuations of the refractive index. The wave front of radiation traveling through the atmosphere changes as it encounters inhomogeneities in the refractive index, degrading optical image quality. The intensity of the turbulent fluctuations of the atmospheric refractive index is described by the refractive index structure function, *C*^{2}_{n}. The maximum telescope resolution is defined by a parameter called *seeing*, which is the integral of *C*^{2}_{n} over the propagation path. Seeing has the units of arc seconds and describes the angle occupied by the star image, at the full-width and half-maximum of its intensity profile (Coulman 1985). The *C*^{2}_{n} and seeing are commonly used by astronomers to describe the turbulent state of the atmosphere at the time of their observations. Roddier (1981) and Coulman (1985) provide reviews of the optical and meteorological aspects of atmospheric turbulence, and the possibility of forecasting atmospheric turbulence from meteorological data was raised by Coulman et al. (1986).

Adaptive optics (AO) allows for a partial correction of the image degradation caused by atmospheric turbulence. See Beckers (1993) for a review. There are a number of physical limitations to AO performance, and successive generations of more sophisticated techniques have been developed. A numerical method that predicts *C*^{2}_{n} and seeing can be used by astronomers to help optimize AO operation. In addition, numerical forecasts of atmospheric turbulence can serve as guidance to schedule telescope usage, since different turbulence conditions suit different types of observations (Bougeault et al. 1995; Masciadri 1998; Masciadri et al. 1999a, b; Businger et al. 2002). Numerical forecast model predictions of the vertical distribution of atmospheric turbulence are also useful for site monitoring by allowing the climatology of atmospheric turbulence and seeing to be constructed.

The first attempt to simulate atmospheric turbulence for astronomical applications using a mesoscale model was made by Bougeault et al. (1995). A turbulence production module was incorporated into the hydrostatic Météo-France limited-area numerical weather prediction model, Péridot (Imbard et al. 1986). The results were promising but some limitations were found, mostly due to the nonhydrostatic assumption and low horizontal resolution of the model (Bougeault 1997). Masciadri (1998) and Masciadri et al. (1999a, b) were able to successfully produce three-dimensional *C*^{2}_{n} maps and integrated astroclimatic parameters, using a nonhydrostatic model with a horizontal resolution of 500 m. Subkilometer resolution proved to be an important factor in the reconstruction of better optical turbulence maps (Masciadri 1998). Validation and calibration of the optical turbulence algorithm implemented in the Meso-NH have led to substantial improvements (Masciadri and Jabouille 2001; Masciadri et al. 2004).

In this paper, we describe an atmospheric turbulence prognostic algorithm implemented in the fifth-generation Pennsylvania State University–National Center for Atmospheric Research Mesoscale Model (MM5), which is currently run operationally at the Mauna Kea Weather Center (MKWC). Output from the algorithm provides guidance for MKWC forecasters and the Mauna Kea astronomy community.

Although there are similarities between the approach described here and that of Masciadri (1998) and Masciadri et al. (1999a, b, 2004), there are also significant differences. In particular, a different turbulent scheme is used and the horizontal and vertical model resolution used here differs from those in the Meso-NH model used by Masciadri. Our initial objective was to implement the optical turbulence algorithm without compromising the current operational schedule; in other words without changing the current model grid resolution.

This paper describes a first step toward designing and implementing an algorithm to estimate future average *C*^{2}_{n} and seeing in the context of an operational mesoscale numerical weather prediction (NWP) model. The algorithm must address the problem that optical turbulence is a subgrid-scale process for NWP models. The goal of this research is to develop and implement a *C*^{2}_{n} parameterization scheme that provides an estimate of the average contribution of optical turbulence to each model grid point from the surface to the top of the model atmosphere. The method does not provide a prognosis of the high-frequency variability that characterizes optical turbulence observations; rather it calculates an average *C*^{2}_{n} estimate over the model grid scale. In the future, the algorithm will be calibrated using summit observations and refined as numerical model grid resolution improves.

## 2. MKWC forecast system overview

The MM5 used at the MKWC is a nonhydrostatic, primitive equation model with a terrain-following coordinate (Grell et al. 1985). It has multiple nesting capabilities to enhance the simulation over the area of interest. The operational model configuration encompasses four two-way nested domains, with horizontal grid spacing of 27, 9, 3, and 1 km, corresponding to a time steps of 81, 27, 9 and 3 s for the four domains, respectively (see Fig. 1). The inner domain is centered on Mauna Kea’s summit, which has an elevation of 4215 m above sea level. Terrain data with a resolution of 30″ (∼0.9 km) obtained from the U.S. Geological Survey are used.

Forty levels in the vertical are used. The height above the ground of the first vertical level is 14 m. The vertical spacing is on the order of tens of meters for the levels nearest the ground, where the greatest contribution to seeing occurs. The spacing gradually increases with height as shown in Fig. 1. The MM5 model uses a staggered vertical-levels system that allows more accurate calculation of the vertical derivatives than when a staggered system is not employed. Given the analytic form of the *C*^{2}_{n} algorithm and the way the integration is performed, *C*^{2}_{n} on the two highest levels is not calculated. Therefore, *C*^{2}_{n} is calculated up to ∼15 km.

The vertical structure used in this study differs from the one used in the work of Masciadri et al. (1999a) and Masciadri et al. (2004). In their work the first grid point is 50 m above the ground; a stretching of 30% in the spacing is applied within the first 3 km, above which the resolution is held constant at 600 m to a height of 20 km. Water vapor and temperature fluctuations tend to be larger in the lowest atmosphere, so the contribution to seeing above 10 km is generally quite small. However, high-altitude turbulence caused by free atmosphere shear (often associated with jet streams) or induced by the underlying complex terrain, can significantly affect the wave front phase de-correlation angle. Therefore, high-altitude turbulence also affects the corrected field of view of current-generation AO systems (Beckers 1993) and numerical prediction will benefit from increases in model vertical grid resolution aloft. A sensitivity study regarding the impact of changes in vertical resolution on seeing is the subject of future work.

The MM5 physics package used in this study includes (i) the grid-resolvable Reisner-2 moisture scheme that models graupel and ice condensation nuclei and allows coexistence of mixed water phases (Reisner et al. 1998), (ii) the Kain–Fritsch cumulus convection scheme (Kain and Fritsch 1990) for the 27-, 9-, and 3-km resolution domains, (iii) a Gayno–Seaman turbulence kinetic energy (TKE)-based planetary boundary layer (PBL) scheme (Gayno 1994; Shafran et al. 2000), and (iv) a longwave–shortwave radiation scheme that allows interaction with water vapor, clouds, precipitation, and the surface (Stephens 1978; Garand 1983).

Local high-resolution analyses of the atmosphere over the North Pacific area, created by the Local Analysis and Prediction System (LAPS), are used to initialize MM5 (McGinley 1989; Cherubini et al. 2006). Boundary conditions are updated every 6 h using model output from the National Centers for Environmental Prediction Global Forecasting System.

## 3. The physics of *seeing*

The turbulent fluctuations of the atmospheric refractive index *n* along the direction *r* are described by the refractive index structure *D _{n}*(

*r*):

Angle brackets in (1) indicate the ensemble average. For locally isotropic turbulence fields,^{1 }Kolmogorov (1941) has shown that

where *C*^{2}_{n} is called the refractive index structure coefficient, which can be considered a measure of the strength of turbulence. An extensive review of this formulation can be found in Roddier (1981).

Since the atmospheric refractive index *n* is a function of the air temperature *T* and concentration of water vapor *q*, the fluctuation of the refractive index is a function of the temperature fluctuation, water vapor fluctuation, and the fluctuation of the combined temperature and humidity fields. Fluctuation of water vapor is negligible at optical wavelengths while observations from Friehe et al. (1975), Antonia et al. (1978), Weseley and Hicks (1978), and Coulman (1980) have shown that the combined effects of humidity and temperature can be significant in a marine boundary layer environment. Because the marine boundary layer resides below Mauna Kea’s summit most of the time, calculation of *C*^{2}_{n}, can be reduced to the calculation of the atmospheric temperature structure coefficient, here expressed in terms of potential temperature *θ*, by means of the following relation, which is valid under the assumption of pressure equilibrium (Bean and Dutton 1966; Masciadri and Jabouille 2001):

An important parameter in optical astronomy is the Fried number or coherence length *r*_{0}, obtained by vertical integration of the refractive index structure coefficient *C*^{2}_{n} along the optical path *L*:

where *λ* is the wavelength at which the telescope is observing. The Fried parameter represents the equivalent diameter of a telescope whose angular resolution is not strongly limited by the atmospheric optical turbulence.

The image quality of a telescope is then described by an optical parameter called the *seeing* of the telescope, *ϕ*. Seeing is defined as the full width at half maximum of a star image at the focus of a large diameter telescope and, for an observation at zenith, seeing is related to *r*_{0} (and therefore *C*^{2}_{n}) by

In this paper *λ* = 0.5 *μ*m, chosen as a representative wavelength for optical astronomy.

Following Coulman et al. (1986) the temperature structure coefficient is calculated according to

where ɛ* _{θ}* is the molecular destruction rate of temperature variance for which we use the parameterization formulated in Bougeault et al. (1995), which follows the study of André et al. (1978),

and ɛ is the energy dissipation rate of the kinetic energy in the velocity field. Equation (6) assumes that the power spectrum of *θ* follows Kolmogorov’s law (Kolmogorov 1941), which is valid in the inertial subrange. Equation (7) assumes that the contributions from the triple correlation and the radiative dissipation in the turbulent energy budget equation can be neglected (André et al. 1978). Therefore, the problem of calculating *C*^{2}_{n} is reduced to calculating the temperature turbulent fluxes and the energy dissipation rate, ɛ. Both these quantities are calculated in the boundary layer parameterization scheme employed in our configuration of the MM5 model. Despite the name, the boundary layer parameterization scheme solves the turbulent kinetic energy budget equation and calculates the temperature turbulent fluxes throughout the entire atmosphere not just within the boundary layer. Expressions for both and ɛ are given in the following subsection, where the TKE-based boundary layer scheme is presented.

The reader will find that the physical parameterization used for the energy dissipation rate, ɛ, and the temperature turbulent fluxes, , are somewhat different from those used in Masciadri et al. (1999a). By using the values for and ɛ calculated by the model, consistency between the seeing algorithm and the rest of the model is maintained.

### a. Introduction to the Gayno–Seaman scheme

The MM5 model offers a wide choice of PBL parameterization schemes (see http://www.mmm.ucar.edu/mm5 for more details). They can be divided into two classes, the “profile based” schemes and the TKE-based schemes. The Gayno–Seaman scheme is a TKE-based parameterization scheme built on work by Mellor and Yamada (1982). This scheme solves the prognostic equation for TKE, or *E*, in which the dissipation of TKE into heat by molecular viscosity is defined as

where *τ*_{0} is the dissipation rate defined as (Gayno 1994; Ballard et al. 1991)

In (9) *c*_{0} and *F*_{L} are dimensionless constants and *F*_{L}*l*^{2}*N* ^{2} is an asymptotic stability term, with *N* being the buoyancy frequency. The basic length scale *l* is given by the Blackadar formula

where *k* is the von Kármán constant and *λ*_{m} is a mixing length defined as one-third of the depth of the turbulence layer containing the height *z*. The depth of the turbulence layer is calculated as the layer where the Richardson number is less than the critical Richardson number.

The temperature turbulent fluxes in the Gayno–Seaman scheme are described in terms of eddy diffusivity,

where *γ*_{g} is a countergradient heat flux term added in the scheme to correct the vertical eddy flux in convective conditions in the boundary layer only (see appendix A) and *K*_{h} is the eddy diffusivity function of TKE, calculated according to

where *l*_{h} is the mixing-length scale that can be expressed as a nonlinear function of the basic length scale *l*, TKE, the buoyancy term *N* ^{2}, and the velocity shear *S*^{2} (see appendix B).

## 4. Instrumentation and optical turbulence data

The vertical distribution of turbulence over Mauna Kea was measured as a part of a site characterization campaign held during October and December 2002. Generalized Scintillation Detection and Ranging (G-SCIDAR) and Multi-Aperture Scintillation Sensor (MASS) instruments operated jointly from 21 to 24 October. Only G-SCIDAR was operated for the December portion of the campaign, which extended from 12 through 18 December.

G-SCIDAR is an instrument that remotely measures the vertical distribution of atmospheric turbulence (Vernin and Roddier 1973; Avila et al. 1997). G-SCIDAR analyzes stellar scintillation of a binary star target produced by turbulent layers present in the atmosphere. This technique relies on the user being able to locate a suitable binary pair in the direction of interest and requires an aperture size of the feeding optics greater than 1 m. For this reason, a G-SCIDAR is not an ideal instrument for regular/operational optical turbulence monitoring. During the 2002, Mauna Kea campaign a G-SCIDAR was installed at the University of Hawaii (UH) 88-in. telescope. Turbulent profiles were determined with 20-s integration time. The vertical resolution of the *C*^{2}_{n} profiles measured by the G-SCIDAR depends on the double stars’ separation (Vernin and Azouit 1983).

Tokovinin et al. (2005, their Table 1) lists the double stars observed during the Mauna Kea campaign and the corresponding vertical resolutions that are of the order of hundreds of meters in the lowest atmospheric layers and thousand of meters at ∼20 km above the ground. The G-SCIDAR instrumental noise at any particular level is equivalent to ∼10^{−18} m^{−2/3}. This leads to an error in the Fried parameter *r*_{0} and therefore in seeing of ∼5% (Masciadri et al. 2002; Prieur et al. 2001, 2004). G-SCIDAR measurements during the 2002 Mauna Kea campaign show that half of the turbulence integral is caused by turbulence below 0.7 km and that this contribution is uncorrelated with the free-atmosphere turbulence (Tokovinin et al. 2005).

On the other hand, MASS relies on the analysis of scintillation of single stars (Kornilov et al. 2003) and provides an integrated value of *C*^{2}_{n} over six atmospheric layers. These layers are centered at 0.5-, 1-, 2-, 4-, 8-, and 16-km altitude above the site. MASS response functions are close to triangular (Fig. 2), so, for example, the turbulent intensity in the 2-km MASS layer is the integral of the optical turbulence profile multiplied by the response function that is zero at 1 km, reaches one at 2 km, and drops back to zero at 4 km.

During the Mauna Kea campaign MASS was installed at the 24-in. UH telescope. The telescope was pointed to a bright single star near zenith and a 1–2-h sequence of 1-min measurements was started. MASS is insensitive to turbulence near the ground, for example, below 250 m (Tokovinin et al. 2005). It can only measure seeing in the atmosphere above the surface layer. The relatively small size of the instrument makes it versatile for regular observation; however, it offers only a low-resolution *C*^{2}_{n} profile. The MASS instrumental error in *C*^{2}_{n} is estimated to be ∼5% of the integrated seeing value and ∼10% for any one layer (Tokovinin et al. 2005).

In the literature, the term *seeing* usually refers to the total seeing or the integral of *C*^{2}_{n} from the ground through the entire atmosphere. In this paper, we will hereinafter refer to seeing as the free seeing measured by MASS, which ignores the surface layer contribution to seeing. All observations used for comparison in this paper have been corrected for zenith angle.

Quantitative comparisons of turbulence profiles measured by G-SCIDAR and MASS during this campaign are presented and discussed in Tokovinin et al. (2005). A good correlation between the two independently calibrated instruments is found in the highest levels. Larger discrepancies and errors are documented for MASS measurements in the lower layers, particularly when multiple layers of comparable intensity turbulence are present. MASS reliability is greater when one dominant turbulent layer is present.

To estimate this systematic error and compare the measurements from the two instruments, the G-SCIDAR seeing was calculated by integrating the G-SCIDAR *C*^{2}_{n} over the six MASS layers, and then the MASS seeing and the G-SCIDAR seeing data were averaged over 5-min intervals. The comparison results in an RMS of 0.1 arcs for the nights of 21–24 October, when both instruments were operating. The average absolute value for seeing on these nights ranges from 0.4 to 0.6 arc s. Therefore, the systematic relative error is from ∼15% to 20%. In a level-by-level comparison of the two instruments, larger systematic errors occur at some levels in the MASS measurements.

The angles to zenith during the G-SCIDAR observations were large because of the lack of suitable stars. The G-SCIDAR observed objects within 60° of zenith. The observed objects for the MASS were generally different from the G-SCIDAR and were limited to zenith distances of less than 45°. During part of the nights of 21–24 October, MASS and G-SCIDAR looked at stars in the same constellations. Specifically, the MASS looked at *α*-Aries and G-Scidar looked at *γ*-Aries, which are separated by 5.2°, between 1 and 3 h every night. The measurements taken during the part of the nights during which the two instruments observed in the same direction within the small 5.2° field of view have been selected to calculate the systematic relative error in observations of *C*^{2}_{n} between the two instruments on each MASS level. This approach minimizes the impact of any horizontal nonuniformity of the turbulence on the measurements. The relative errors between the two instruments are quite large in the lowest three MASS layers, on the order of 50%–100%. Smaller systematic errors, on the order of 30%–70%, are found in the three highest levels.

These findings have to be kept in mind when comparing the model results to available measurements. At this time, our ability to determine the model accuracy is limited by the systematic errors found among the observations.

## 5. Case studies and results

Once the algorithm to calculate *C*^{2}_{n} and seeing was implemented in the MM5 code, the model was rerun for several nights during the 2002 Mauna Kea monitoring site campaign, including the nights of 22 and 23 October and 15 and 16 December.

For each of the nights under investigation, MM5 was initialized at 0000 UTC of that same nominal day. Turbulence data from both MASS and G-SCIDAR were collected between 0600 and 1600 UTC every night. Therefore, the forecast output used in the comparison is from the 6th to the 16th hour of MM5 forecast. The *C*^{2}_{n} simulated data are from the MM5 innermost domain, with horizontal resolution of 1 km. Although the implemented algorithm produces *C*^{2}_{n} in the surface layer, only simulated data from 250 m and up have been used in the comparison to match the MASS observation range.

The MM5 optical turbulence algorithm calculates *C*^{2}_{n} at each model level and for every model time step. The model time step used for the innermost domain is 3 s. Using the MASS triangular response functions *wf* (*z*) [Eq. (13) and Fig. 2], turbulence integrals for six atmospheric layers are calculated for G-SCIDAR and MM5-derived data:

where the layer heights *h*_{MASS}= 0.5, 1, 2, 4, 8, and 16 km (after Tokovinin et al. 2005).

Turbulence integral calculations allow for an easy comparison of the model-simulated turbulence conditions with MASS and G-SCIDAR observations. Time series of observed and simulated seeing and average *C*^{2}_{n} profiles at the summit will also be compared below.

The weather pattern for 22 and 23 October 2002 was dominated by a zonal upper-level flow and a surface ridge located north of the state of Hawaii (Fig. 3). A weak midlevel ridge was strengthening, following passage of an upper-level trough that swept through the central Pacific area a few days earlier. Consequently, the trade wind inversion was gradually strengthening, winds aloft were moderately strong, and substantial shear was present between the summit and the upper troposphere (Fig. 4). At the summit, precipitable water ranged between 3 and 4 mm on 22 October and between 2 and 3 mm on the following night. Summit wind speed was in the range from 2.5 to 5 m s^{−1} on both nights. Under these conditions, there are contributions to Mauna Kea summit seeing from turbulence at the upper levels related to high wind shear and moderate contributions from friction at the lowest atmospheric levels. Only one synoptic time is shown in Fig. 3, because the synoptic conditions for 22 October were relatively persistent.

The MM5 forecasts for both the October and the December case were quite good in that the model forecast errors of summit variables and large-scale pattern were well within the model RMS (Cherubini et al. 2006). Therefore, discrepancies between predicted and observed seeing are not the result of a poor model weather forecast for the nights under investigation.

By plotting the *C*^{2}_{n} at several levels, one can discern which atmospheric level is contributing most to the seeing (Figs. 5 and 6). Given the atmospheric conditions present, optical turbulence was consistently observed within the boundary layer and at the highest levels. A comparison between the signals recorded from the G-SCIDAR and MASS illustrates the difficulty in obtaining an absolute measure of *C*^{2}_{n}. A reasonable degree of agreement between G-SCIDAR and MASS data is found for the highest levels while a certain amount of overestimation by MASS is evident at the 0.5-km level. The discrepancies found between the MASS and G-SCIDAR datasets are explained in detail in Tokovinin et al. (2005).

As expected, when MM5 output is compared with the observed data for the nights in October, the model is not able to reproduce the high-frequency variability shown by the observed data. The model output becomes even smoother at higher altitude. This behavior is understandable in terms of model resolution. A horizontal grid spacing of 1 km is coarse relative to the scale of turbulence. The vertical spacing, although reasonably fine in the lowest layers of the troposphere, is of order 600 m at ∼8 km above the summit, and increases even further near the top of the model.

Following Tokovinin et al. (2005), given the magnitude of the variability in the observed data particularly for the MASS dataset and the large temporal variability of optical turbulence, it makes sense to compare averaged quantities (Table 1). The average seeing for each night at each level was calculated using the 1-h average seeing values from 0600 to 1600 UTC. Each 1-h average value was included in the nightly average only if more than 20 profiles were available from the instruments. This latter condition was imposed to ensure confidence in the hourly averages, since the observations naturally present gaps through the night.

Table 1 shows that the MM5 model output for both nights overestimates the turbulence for the levels centered at 1 and 2 km above the site. Nevertheless, the model predicts average values of seeing that are quite close to those observed on both nights (Figs. 5 and 6).

Average profiles of *C*^{2}_{n} from the G-SCIDAR for 22 and 23 October 2002 show two turbulent layers centered at about 4 and 12 km above the site and a boundary layer in which *C*^{2}_{n} rapidly increases nearer the ground (Fig. 7). MM5-derived *C*^{2}_{n} average profiles also show two levels of upper-level turbulence, though they are slightly displaced from the ones observed. Both profiles show a significant contribution to the optical turbulence from the lowest 3 to 4 km on 22 October 2002.

The following approach was used for all nights analyzed to construct the average G-SCIDAR profile in an effort to retain some of the higher-resolution character of the individual profiles. The Fried number, *r*_{0}, for each individual profile is calculated first; a subset of the profiles, consisting of 10% of the total number in which *r*_{0} is nearest the median r_{0}, is then selected. The average profile is calculated using this subset only. This approach provides a better representation of the individual profiles than a simple layer-averaged profile approach, which washes out sharp features.

Knowledge of the nightly averaged profile of predicted *C*^{2}_{n} can help astronomers in calibrating their AO systems. The calibration requires estimating a priori where most of turbulence is expected to occur.

A well-developed ridge was centered over the state of Hawaii on 15 and 16 December (Fig. 8). The ridge resulted in strong subsidence, a well-defined trade wind inversion (Fig. 9), relatively warm temperatures, and low precipitable water (∼1 mm) at the summit level. Winds aloft were light and shear between the summit level and the upper levels of the atmosphere was weak. Summit wind speeds ranged between 5 and 10 m s^{−1} on both nights. Under these conditions, smaller contributions to optical turbulence are expected from the highest atmospheric levels than from the lowest levels.

During the nights of 15 and 16 December, there is reasonable agreement between the observed and simulated *C*^{2}_{n} in the lower troposphere (0.5, 1, 2, 4 km) (Figs. 10 and 11). A certain degree of overestimation is evident at the two highest levels, particularly for the night of 16 December (Fig. 11). In spite of the model overestimation of *C*^{2}_{n} at the upper levels, predicted seeing (total) values are reasonably well correlated with the observed seeing (Table 2).

It is suggested that high temporal variations in simulated values of *C*^{2}_{n} at the upper levels on the night of 16 December may be explained by the poor vertical resolution at higher altitudes and the model’s known tendency, under certain conditions, to show numerical instabilities at the top model boundary as a result of the interaction of the atmospheric flow with complex orography (Klemp and Durran 1983; Gallus and Klemp 2000; Klemp et al. 2003). Future work will investigate possible solutions to this problem.

G-SCIDAR profiles show weak turbulence above 10-km altitude and a rapid but uniform increase in *C*^{2}_{n} with decreasing altitude (Fig. 12). The MM5-simulated average profile of *C*^{2}_{n} correlates reasonably well with the observed profile up to 5-km altitude. A discrepancy between the two profiles is evident at higher altitudes (above 6 km). Nevertheless, the agreement between observed and predicted seeing for both nights is reasonable, consistent with the fact that the lowest tropospheric levels contribute the greatest amount to the seeing.

It is clear that a good correlation between measured and predicted seeing can result from a less good correlation between measured and predicted *C*^{2}_{n} profiles, as seems to be the case for the December nights under investigation. A similar problem is also found when comparing *C*^{2}_{n} profiles from different instruments (i.e., MASS versus G-SCIDAR). Nevertheless, predicting seeing is an important component in predicting the quality of the observing conditions at the summit of Mauna Kea. Therefore, there is considerable motivation to further understand and reduce the discrepancies between measured and predicted profiles through this and future research.

Output from MM5’s optical turbulence algorithm also provides the opportunity to construct three-dimensional forecasts of seeing and *C*^{2}_{n} (Fig. 13). The October case illustrates an elevated turbulence layer at ∼300 hPa (Fig. 13a), whereas turbulence in the December case is primarily distributed nearer the surface (Fig. 13c). The horizontal maps of seeing show increasing seeing values on the lower slopes of Mauna Kea, associated with clouds and background turbulence in the planetary boundary layer beneath the trade wind inversion (Figs. 13b,c). Output like that shown in Fig. 13 provides useful guidance for forecasters, and archived model data are potentially useful for site characterization efforts.

## 6. Conclusions and future work

In this paper, we describe an atmospheric turbulence algorithm that has been implemented in the mesoscale weather prediction model MM5, which is currently run operationally to provide short-range weather forecasting guidance at the Mauna Kea Weather Center (MKWC). The horizontal grid spacing of MM5 over the summit of Mauna Kea is fixed at 1 km, whereas the vertical resolution is density weighted, coarser close to the ground (∼tenths of meters) and increasing to order ∼1 km from the tropopause up to the top of the model, fixed at 10 hPa (∼25 km). Given the current model configuration, optical turbulence must be regarded as a subgrid-scale process. The optical turbulence algorithm provides an estimate of the turbulence at the model grid point from subgrid-scale optical turbulence processes. This scheme makes use of information regarding turbulence kinetic energy (TKE) and temperature turbulent fluxes parameterized in the model to calculate the turbulent fluctuations of the atmospheric refractive index, *C*^{2}_{n}, and *seeing*.

The vertical distribution of turbulence over Mauna Kea was measured as part of a site characterization campaign during several nights in October and December 2002. The MM5, modified with the optical turbulence algorithm, was rerun for 22 and 23 October and for 15 and 16 December.

Overall, the implemented algorithm was able to reproduce the average observed behavior of the observed *C*^{2}_{n}, particularly for the lowest atmospheric level. As a result, for the studied cases the predicted and observed average values of seeing are reasonably well correlated. However, the observed data show a temporal variability that the model is not able to capture, particularly at the highest levels where model vertical resolution degrades.

Operationally, each MM5 model run takes 4 h and 40 min on average for a 60-h forecast. The effect of the implementation of the *C*^{2}_{n} algorithm within the model on the simulation length is negligible because the *C*^{2}_{n} algorithm makes use of variables that have already been calculated by the model for other purposes. The seeing and *C*^{2}_{n} plots available on the MKWC Web page are produced in postprocessing and they are available at the same time as plots of other meteorological variables.

### a. Future work

The purpose of this paper is to present the details of the seeing algorithm, using specific cases to illustrate our approach. It is clear that a sample of four cases is insufficient to undertake a statistical analysis of the success of the method. For a robust validation of the model, many additional cases are needed. In Cherubini et al. (2007, manuscript submitted to *J. Appl. Meteor. Climatol.*), seeing observations routinely provided by several telescopes on Mauna Kea will be used to validate and calibrate the seeing algorithm following a statistical approach. Good versus poor seeing forecast periods will be investigated relative to their corresponding weather conditions and to the overall accuracy of the MM5 weather forecasts. A poor seeing prediction might in fact be a consequence of a poor model weather forecast rather than a deficiency in the seeing algorithm performance. The problem of predicting *C*^{2}_{n} in the lowest model levels will also be addressed in this second paper and the *C*^{2}_{n} algorithm will be tuned by calibration of the eddy diffusivity in the atmospheric layers closest to the ground. Daily forecasts posted and archived at our Web site provide a preview of the results that will be included in this follow-on paper (see http://mkwc.ifa.hawaii.edu/).

In future studies (i) data from slope detection and ranging (SLODAR), an instrument able to accurately measure optical turbulence in the lowest atmospheric layers, will be analyzed to improve *C*^{2}_{n} parameterization in the lowest atmospheric levels, (ii) the sensitivity of the algorithm performance to the horizontal and vertical resolution of the model will be undertaken, and (iii) numerical instabilities associated the model atmospheric flow interaction with complex topography will be investigated.

## Acknowledgments

We thank Aziz Ziad and Jean Vernin (Université de Nice) for the SCIDAR instrument data and Andrei Tokovinin (Cerro Tololo Interamerican Observatory) for the MASS instrument data. We are grateful to Dr. E. Masciadri for her thoughtful reviews of this paper in its various stages. We also thank two anonymous reviewers for their comments. The National Astronomical Observatory of Japan provides computer resources, which are located at the Subaru Telescope facility in Hilo, Hawaii, for the MM5 simulations.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

_{n}

^{2}simulations.

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

### APPENDIX A

#### The Countergradient Heat Flux Term γg

The analytic form of the countergradient heat flux term is given by the following relationship (Gayno 1994; Shafran et al. 2000):

valid for *z* < 1.2 *h*_{PBL}, where *h*_{PBL} is the local height of the planetary boundary layer.

In (A1) *H _{s}* is the surface heat flux, while

*w** is the convective vertical velocity scale given by

In (A2) *θ _{υ,}*

_{0}is the virtual potential temperature in the lowest model level and

*g*is the acceleration of gravity.

### APPENDIX B

#### Length Scales Used in Coefficients of Turbulent Transport

The length scales *l _{θ}* and

*l*in the diffusion coefficients

_{m}*K*and

_{h}*K*in section 3 are given by the following expressions:

_{m}and

where

Here, *c _{i}* are empirically determined nondimensional constants relating all dissipation and return-to-isotropy time scales in the diagnostic algebraic expressions for the second-order equation for the turbulent kinetic energy so that

*τ*[the relation

_{i}= c_{i}τ*τ*

_{0}= c

_{0}

*τ*is given in section 2 of Ballard et al. (1991)]. Here,

*c*

_{0}= 5.524 27,

*c*

_{1}= 1.046 94,

*c*

_{2}= 1.413 37,

*c*

_{3}= 3.084 92, and

*F*= 7.097 59. From Launder (1975),

_{L}*a*

_{2}= ⅓. (Adapted from Ballard et al. 1991.)

## Footnotes

*Corresponding author address:* Dr. Tiziana Cherubini, Department of Meteorology, University of Hawaii at Manoa, 2525 Correa Rd., Honolulu, HI 96825. Email: tiziana@hawaii.edu

* School of Ocean and Earth Science and Technology Contribution Number 7258.

^{1}

This formulation is valid when *l* ≪ *r* ≪ *L,* with *l* being the inner scale, or the size below which viscous effects are important and energy is dissipated into heat, and *L* being the outer scale or the size above which isotropic behavior is violated.