3.1

Spatial and temporal resolutions

Different kinds of NWP software for modelling the atmosphere are available, such as general circulation models (GCM), mesoscale models or large-eddy simulations (LES).²⁵ Each of them resolves atmospheric phenomena at different scales, offering their own spatial and temporal resolutions.

Since WRF belongs to the category of mesoscale models, it is applied to a limited area of Earth and can simulate mesoscale atmospheric phenomena ranging from kilometers to hundreds of kilometers. As mentioned previously, grid nesting is used to increase the spatial resolution near the location of interest, achieving horizontal resolution of 1 km.²⁰ This limitation is mostly computational and the resolution of the orographic data must also be considered when reducing further the grid spacing.

Vertically, the grid spacing is usually varying non-linearly, with smaller distance between levels close to ground and larger distance between levels at the top of the atmosphere. Close to ground, the spacing can be on the order of 10 meters, whereas, at high altitude, the spacing can be more than 500 m. In order to further reduce the vertical grid spacing, more levels should be added in the NWP software at the price of a larger computation time. In Ref. 20, it is suggested to use grid spacing of 100 meters or less to resolve OT.

Regarding the temporal resolution, it is intrinsically related to the time step required for avoiding numerical instabilities in NWP simulations (e.g. 1 to 5 seconds²⁰). However, this is not the temporal resolution at which output meteorological data are recorded. Instead, they are saved every 5 minutes in this work.

Following these considerations, the configuration for WRF simulations at Redu is given in Tab. 1. The three nested domains are named d01, d02 and d03, with d01 being the largest but coarsest domain and d03 being the smallest but finest domain. In each domain, 250 vertical levels are used leading to an average vertical spacing between levels of approximately 80 meters, from the ground to 20 km of altitude. Levels are almost equally-distributed, with the largest level spacing being 100 meters and the smallest one being 60 meters. The top pressure of the atmosphere set for WRF simulations is 5000 Pa. For comparison with measurements, simulations start at least 6 hours prior to the time window where the measurements are available (6 hours of lead time).

Table 1:

Grid parameters from WRF simulations.

The configuration of the WRF physical schemes used is: WSM6²⁶ for microphysics, Tiedtke²⁷ for cumulus physics, Dudhia²⁸ and RRTM²⁹ for shortwave and longwave radiations, and revised MM5³⁰ for the surface layer. For the planetary boundary layer physics, the MYNN 2.5³¹ scheme is used as it solves for the turbulent kinetic energy (TKE) that can be added to WRF simulation outputs. This TKE is then used in the chosen model.

3.2

Choice of model

The choice of model is also a challenge. It corresponds to the step that enables to go from meteorological data to profiles useful to get OT-related quantities. The literature about optical models is rich and different models have been developed throughout the years. Several classifications of these models are possible:

where a²α are constant parameters, L₀ is the outer scale of turbulence and is the refractive index vertical gradient. Models for L₀ and M are then required.

Alternatively, a last category of models has recently emerged, grouping numerical models solving for the TKE in fluid mechanics simulations. They are of particular interest for OT forecast as meteorological data are often obtained thanks to NWP simulations. Hence, this latter can be adapted to compute and extract the TKE and then use it in a model.^{18, 20, 36}

Furthermore, other classifications are also possible. For example, models can also be classified depending on the part of the atmosphere they describe (free atmosphere or surface layer), on the location they have been designed for (e.g. for observatories at high altitude), on the time of the day (daytime or nighttime models), on the hypotheses used, as well as depending on their primary utility (for astronomical observations or for optical communications) and on the targeted wavelengths. Therefore, all of these factors must always be taken into account when using any model.

As a summary, the choice of model for OT forecast is not trivial. For OT characterization and prediction based on NWPs, chosen models are mostly parametric such that they depend on local meteorological quantities. This dependence can either be expressed through analytical expressions (theoretical models) or through expressions involving the TKE (numerical models). They usually describe the whole atmosphere, even though they can use a separate model for the boundary layer. They can also include calibration parameters based on empirical measurements at the location of interest.

In recent works, numerical models solving for the TKE have been used, e.g. in Meso-Nh^{18, 37} or with WRF.^{20, 38, 39} Theoretical models based on Tatarskii’s equations are also commonly used.^{11, 15, 40} Finally, empirical models are also sometimes considered, such as the Trinquet-Vernin model.^{41, 42}

In this study, the model, Astro-Meso-Nh, presented in Ref. 18 has been chosen, namely because it is parametric and has already been used extensively in OT prediction for astronomy.^{18, 37} It relies on the TKE that is directly extracted from NWP simulations. As for several models, it is based on Gladstone’s relationship³³ linking the refractive index structure parameter to the temperature structure parameter :

where p is the pressure in hectopascal and T is the temperature in kelvin. Equation (2) is obtained for a wavelength of 500 nm and is assumed to be valid for all visible wavelengths.⁷ Different expressions for exist but they are usually derived from the Tatarskii’s expression giving as the product between and the square of the temperature vertical gradient.² In Masciadri,¹⁸ the chosen model is

with the grid-averaged value of the potential temperature.¹³ In stable layers, the parameter ϕ₃ is equal to 0.78 and the mixing length L is the Deardoff length,⁴³ i.e.

where e is the TKE, g is the gravity of Earth and θ_v is the virtual potential temperature⁴⁴ obtained from θ_v = θ (1 + 0.61r) for unsaturated air with mixing ratio r of water vapor.

Hence, this model involves macroscale meteorological quantities (pressure p, temperature T, potential temperature θ and mixing ratio r) that are directly obtained from WRF simulations. Computing the mixing length also requires the turbulent kinetic energy e that is thus added to the outputs of WRF. A similar model is used in Ref. 38, also relying on the Tatarskii’s expression and a TKE model for the mixing length (assumed to be the outer scale L₀ in their work).

Figure 2 illustrates the contribution of the different factors in Eqs. (2) and (3) for Redu at 21:00 UTC on February 14, 2019. Turbulent layers, i.e. large values in the profile of Fig. 2f, are associated with large potential temperature gradient values (Fig. 2c) or large mixing length values (Fig. 2d). As seen in Fig. 2b, turbulent layers are associated with peak values of TKE. The TKE has a minimum value of 5 × 10^–5 m²s^–2 set by default in WRF. A possible calibration of this minimum value based on profile measurements has already been presented in previous work³⁶ but was not applied in this study.

Figure 2:

Illustration of the different factors involved in model for Redu at 21:00 UTC on February 14, 2019.

3.3

Boundary layer effects

The boundary layer is the portion of Earth’s atmosphere close to the ground where effects coming from the Earth’s surface are not negligible. Usually, turbulence is quite strong in this layer and is difficult to model. This is particularly a challenge at optical communication sites that are not located at high altitude and that suffer thus from stronger turbulence due to the boundary layer. A possible solution to ease this challenge is the use of hybrid models separating the free atmosphere from the surface layer.⁴⁵

However, in this study, a single model for describing the whole atmosphere is preferred. It is made possible thanks to the use of the TKE that can also provide predictions in the boundary layer.

3.4

Availability of measurements

In order to validate the prediction capabilities of the model, measurements of profiles are not always available and integrated quantities (e.g. seeing, isoplanatic angle, scintillation index, etc.) are often much simpler to obtain. Nevertheless, as previously stated in Sec. 2, knowledge of profiles enables to obtain any integrated quantity of these profiles thanks to analytical expressions.

In the following section, seeing measurements at Redu will be compared to their predictions from NWP simulations. The seeing, noted ϵ₀, corresponds the full width at half maximum of the point spread function obtained when imaging through turbulence. It should be as small as possible and can be used to identify temporal windows during which the AO system will be most effective.²⁵ It is related to the Fried parameter r₀ by ϵ₀ = 0.98λ/r₀, and is an integrated quantity of the profile that is often dominated by the surface layer.⁴⁶ Its analytical expression in the case of plane waves is given by

where h₀ is the starting altitude (above ground) for the integration of the profile and ξ is the elevation angle measured from zenith, assumed to be 0° in the following. Hence, predictions of seeing come from Eq. (5) where simulated profiles have been integrated. They can then be compared to seeing measurements, obtained in this study using a differential image motion monitor (DIMM).⁴⁷

A DIMM is a standard and widely spread instrument for seeing measurement. It can be easily implemented on a small amateur telescope, and uses typically a mask dividing the telescope pupil in two smaller subapertures. The seeing, or Fried parameter, is then estimated from the variance of the differential image motion of the two resulting spots (assuming one observes a star).⁴⁷ In the framework of SALTO,¹⁷ a DIMM instrument has been set up using a Celestron C14 (D = 35.6 cm) with a mask of aperture size of d =12 cm and a baseline of B = 25 cm. The DIMM software has been implemented following the prescriptions of Ref. 47. An illustration of the setup at Redu is given in Fig. 3.

Figure 3:

Illustration of the DIMM used at Redu in the framework of SALTO.

4.1

Results

Currently, three nights of measurements from the SALTO project are available (February 05, February 14 and September 04, 2019) and have been exploited for model validation.¹⁷ Figure 4 depicts the evolution of measured and predicted seeing for February 14, 2019. DIMM measurements are available every 5 to 10 seconds and are relatively noisy (gray curve in the background of Fig. 4). A possible approach to deal with noisy measurements is filtering, this is the reason why a moving-average filter with a window of 50 points has been applied, leading to the orange curve. Alternatively, binning of the measurements in bins of 5 minutes has been performed. The average of all measurements falling in a given bin corresponds to the blue points, and the bars represent their standard deviation. Seeing predictions coming from WRF simulations with the TKE-based model correspond to the green squares, and are available every 5 minutes. Therefore, they can directly be compared with the average seeing value from the binned measurements. Finally, as a reference, seeing predictions coming from integration of the Hufnagel-Valley (HV) 5/7 profile are depicted in red.

Figure 4:

Prediction of seeing on February 14, 2019, at Redu (Belgium) and comparison with measurements.

WRF simulations have been performed with a lead time of 6 hours prior to the first available measurements. Computations of seeing from simulated profiles have been achieved using Eq. (5), starting the integration at an altitude h₀ of 100 meters up to 20 km. This starting altitude approximately corresponds to the first vertical level in WRF simulations and can be seen as a calibration parameter that is further discussed below.

As seen in Fig. 4, realistic seeing values are predicted by the model. However, they do not follow the short-term evolution of the measurements.

Comparison of seeing predictions and measurements for all three nights is depicted in Fig. 5. For every bin of 5 minutes where measurements are available, the measured average seeing is represented on the x-axis, whereas the y-axis shows the predicted seeing value for this 5-minute interval. For perfect seeing prediction, black points should fall on the dashed line. Moreover, the mean and standard deviation of the seeing measured and predicted for a complete night is depicted by the red points and bars.

Figure 5:

Correlation plot between measured seeing and predicted seeing for all three nights at Redu.

Figure 5 shows that even though short-term predictions remain a challenge, long-term trends are similar: when the measured seeing is large for one night, the predicted seeing is also large. There is however an underestimation of the predicted seeing for one night (September 04, 2019), as seen from the cloud of points centered on the red point in (2.6;1). This underestimation results from a particularly low TKE value at the first WRF level, as well as a small gradient of potential temperature at this level.

4.2

Discussion

Inaccuracies of (short-term) seeing predictions have been found to originate from limited prediction capabilities of NWP simulations as well as from the boundary layer .

Prediction capabilities of NWPs Indeed, NWP models solve partial differential equations to obtain meteorological quantities on the desired grid and at given times. Since these equations are quite sensitive to the initial conditions (that are not known perfectly in practice), one can expect discrepancies between predictions and observations. Especially, the time at which a particular weather phenomenon is expected to happen is not always predicted properly, namely for phenomena with a small spatial extension and far in the future. Hence, when considering NWP models to study atmospheric impairments on electromagnetic wave propagation, comparisons are usually performed statistically.⁴⁸

This effect is illustrated in Fig. 6 depicting the predicted seeing and the measurements for all three nights. The predictions come from the TKE-based model applied to NWP simulation outputs, with the only difference that the simulations are performed either with a lead time of 6 hours or a lead time of 12 hours. This variation of lead times introduces some changes in the meteorological quantities predicted during the night, hence impacting the predicted seeing at each instant. Nevertheless, average seeings over a given night remain consistent. A possible solution to reduce such inaccuracies can be the use of mathematical models combining seeing observations and NWP predictions to provide accurate short-term seeing predictions. Such models have already been applied to astronomical sites.¹⁴

Figure 6:

Seeing predictions for all three nights depending on the lead time (6 hours or 12 hours) in WRF simulations.

Boundary layer and seeing For seeing, accurate modelling of the boundary layer is paramount. As seen from Eq. (5), the integration only depends on the profile. Therefore, layers of large values tend to drastically impact the seeing. Since is usually the largest in the boundary layer, this layer dominates the seeing.⁴⁶

As an example, Table 2 provides the average seeing for each night with different upper boundaries for the integral in Eq. (5). In both cases, integration starts from h₀ = 100 m but ends up either at 500 m or at 20 km. Results of Tab. 2 show that seeing at Redu arises mostly from the boundary layer, between h₀ to 500 m. Hence, if one is only interested in seeing predictions, there is no need to know the complete profile. Stated differently, it means that seeing measurements cannot be used solely to validate complete profiles. In order to overcome these limitations, other integrated quantities of profiles can be used. For example, the isoplanatic angle θ₀ involves a weighting function of the altitude in the integration, giving more importance to high-altitude values than to ground values.¹⁵ Multi-aperture scintillation sensor (MASS) measurements can also provide free-atmosphere seeing, above the boundary layer.³⁹ Alternatively, measurements of profiles could be used to validate modelled profiles.

Table 2:

Average seeing and standard deviation per night, integration from h0 to 500 m or 20 km of altitude.

Moreover, extrapolating the value at 100 m (i.e. approximately from the first vertical level in WRF) up to the ground and starting the integration at h₀ = 0 m led to unrealistic seeing values, as large as 6 arcsec. This highlights the need for accurate descriptions of the boundary layer at optical communication sites where OT is largely dominated by this layer. It must be described with a spatial resolution larger than the one used in mesoscale simulations, either relying on hybrid models involving analytical profiles of in the boundary layer,⁴⁵ or using LES to increase the resolution at the price of a higher computational cost.