## OTCliM: generating a near-surface climatology of optical turbulence strength ( C 2 n ) using gradient boosting

Maximilian Pierzyna

Department of Geoscience and Remote Sensing, Delft University of Technology, Delft, The Netherlands

## Sukanta Basu

Atmospheric Sciences Research Center, University at Albany, Albany, USA ∗

## Rudolf Saathof

Faculty of Aerospace Engineering, Delft University of Technology, Delft, The Netherlands (Dated: August 2, 2024)

This study introduces OTCliM (Optical Turbulence Climatology using Machine learning), a novel approach for deriving comprehensive climatologies of atmospheric optical turbulence strength ( C 2 n ) using gradient boosting machines. OTCliM addresses the challenge of efficiently obtaining reliable site-specific C 2 n climatologies, crucial for ground-based astronomy and free-space optical communication. Using gradient boosting machines and global reanalysis data, OTCliM extrapolates one year of measured C 2 n into a multi-year time series. We assess OTCliM's performance using C 2 n data from 17 diverse stations in New York State, evaluating temporal extrapolation capabilities and geographical generalization. Our results demonstrate accurate predictions of four held-out years of C 2 n across various sites, including complex urban environments, outperforming traditional analytical models. Non-urban models also show good geographical generalization compared to urban models, which captured non-general site-specific dependencies. A feature importance analysis confirms the physical consistency of the trained models. It also indicates the potential to uncover new insights into the physical processes governing C 2 n from data. OTCliM's ability to derive reliable C 2 n climatologies from just one year of observations can potentially reduce resources required for future site surveys or enable studies for additional sites with the same resources.

## I. INTRODUCTION

Atmospheric optical turbulence is highly relevant for optical ground-based astronomy and future free-space optical communication (FSOC). Both applications suffer from light getting distorted when propagating through the turbulent atmosphere. In astronomy, turbulent fluctuations of the atmospheric refractive index, known as optical turbulence (OT), cause blurry images and limit the detection of small objects [1]. FSOC links, which use optical beams to transmit data instead of traditional radio waves, experience reduced data rates or even link interruptions due to OT [2, 3]. Therefore, the optical turbulence strength ( C 2 n ) must be carefully considered in the design and operation of such optical systems. This requires robust statistical evaluation of C 2 n over time for the operational sites of interest, known as a C 2 n climatology. Such a climatology, derived from long-term C 2 n data, portrays trends, seasonal variations, or potential anomalies in OT strength.

years at a few carefully selected locations to identify the best one (see, e.g., Hill et al. [4], Sch¨ck o et al. [5]). When envisioning a global FSOC network, site surveys have so far been focused on cloud cover (see, e.g., Fuchs and Moll [6], Poulenard et al. [7], Pham et al. [8]), but similar site surveys targeting OT are expected to become highly relevant in the future.

Obtaining C 2 n climatologies is challenging because turbulence strongly depends on the local topography and varying meteorological conditions. That is why longterm on-site surveys measuring local OT conditions are crucial. For instance, before a new optical observatory is built, OT strength is typically measured for multiple

∗ Also at Department of Environmental and Sustainable Engineering, University at Albany, Albany, USA

However, conducting long-term site surveys at various locations of interest is time-consuming and resourceintensive. While C 2 n can be obtained by post-processing numerical weather model output (e.g., Masciadri et al. [9]), running such models for multiple years is computationally expensive, and the accuracy of the resulting C 2 n is very sensitive to the model configuration [10] and the selected C 2 n parameterization [11]. We address these issues by proposing a novel machine learning-based (ML) approach called OTCliM (Optical Turbulence Climatology using Machine learning). OTCliM aims to extrapolate just one year of measured C 2 n into a multi-year time series, enabling the generation of a comprehensive sitespecific C 2 n climatology with less data than a conventional site survey. As a result, site survey costs can potentially be reduced, climatologies can be obtained faster, and more sites can be surveyed within a given time frame. By leveraging ML, OTCliM can also model complex inputoutput relations, offering an advantage over traditional empirical C 2 n models, which are discussed in the following section.

This study assesses OTCliM's performance using an extensive C 2 n dataset containing measurements from 17

diverse stations. We evaluate the temporal extrapolation capabilities of the trained models and analyze their potential for geographical generalization to other sites. Additionally, we examine the importance of each input variable for predicting C 2 n to probe if the ML models have learned physically plausible relations. That analysis increases confidence in the ML models and could uncover new dependencies of C 2 n from the data.

## II. RELEVANT C 2 n REGRESSION STUDIES

The experimental determination of C 2 n is challenging because it requires expensive instruments, such as scintillometers [12] or meticulous post-processing of high-frequency temperature measurements [13]. Consequently, there has been a sustained effort over several decades to develop parameterizations of C 2 n using more readily available meteorological variables. Conventional approaches to C 2 n parameterization typically involve the application of physics-based analytical equations in conjunction with empirically derived coefficients or regression functions. A popular example is the formulation proposed by Wyngaard et al. [14], based on Monin-Obukhov similarity theory (MOST) [15]. That parameterization called W71 in the following - relates sensible heatflux w θ ′ ′ and friction velocity u ∗ to the strength of temperature fluctuations ( C 2 T ) at height z as

$$C _ { T } ^ { 2 } = \left ( - \overline { w ^ { \prime } \theta ^ { \prime } } / u _ { * } \right ) ^ { 2 } z ^ { - 2 / 3 } g ( \zeta ). \quad \quad ( 1 ) \quad \text{due}$$

Since the refractive index fluctuations quantified by C 2 n are a result of density fluctuations due to temperature fluctuations, C 2 T is linked to C 2 n through the Gladstone equation

$$C _ { n } ^ { 2 } = \left ( A \, P / \overline { T } ^ { 2 } \right ) ^ { 2 } C _ { T } ^ { 2 } \, \text{ \quad \ \ } ( 2 ) \, \text{ \ \ } \text{M}$$

with A ≈ 7 9 . × 10 -5 KhPa -1 for optical wavelengths [16], station pressure P in hPa, and mean air temperature T in K. Function g ζ ( ) in Eq. 1 is a similarity function, empirically determined from observations,

$$g ( \zeta ) & = 4. 9 ( 1 - 6. 1 \zeta ) ^ { - 2 / 3 }, \zeta < 0 \quad \text{(3)} \quad \text{serv} \\ g ( \zeta ) & = 4. 9 ( 1 + 2. 2 \zeta ^ { 2 / 3 } ), \zeta \geq 0 \quad \text{(3)} \quad \text{into}$$

with stability parameter ζ = z/L where L is the Obukhov length L = -u 3 ∗ T/ κg w θ ( ′ ′ ) and z the height above ground. In L g , = 9 81 m s . -2 is the gravitational acceleration of Earth and κ = 0 4 is the Von K´rm´n constant. . a a Various alternative regression-based similarity functions g ζ ( ) were proposed in literature aimed at improving the accuracy and applicability of C 2 n estimates across diverse meteorological and topographic conditions (see, e.g., Savage [17] for an extensive review).

In recent years, different machine learning (ML) techniques have been utilized to derive C 2 n parameterizations.

Similar to W71, ML-based approaches aim to obtain regression models that estimate C 2 n from routine meteorological variables (temperature, pressure, wind speed), gradients (potential temperature gradient or wind shear), or heat and radiation fluxes. However, instead of deriving physics-based functional expressions (cf. Eqns. 1 and 2), which are fitted to observations (cf. Eq. 3), ML directly models the relation between observed C 2 n and the input variables. The power of such ML regression models is that they can model complex multi-variate relations, yielding potentially better C 2 n estimates compared to traditional approaches.

The ML-based C 2 n parameterizations in literature differ in their input variables, the ML regression technique employed, and if C 2 n is measured at a single or multiple vertical levels. Much work has been devoted to singlelevel near-surface C 2 n obtained as point measurements (e.g., from sonic anemometer) or via path-averaging (e.g., from scintillometer). Wang and Basu [18], for example, used fully connected neural networks to estimate singlelevel C 2 n from routine meteorological data and gradients. Vorontsov et al. [19] also utilized neural networks but aimed at estimating path-averaged C 2 n from images of a received laser beam that is broken up into speckles due to turbulence. Jellen et al. [20, 21] employed random forest and gradient boosting machines (GBM) to estimate path-averaged C 2 n in a maritime environment from routine parameters and radiation fluxes. Most MLbased models do not yield analytical equations of the fitted model, unless they explicitly target it as shown by Bolbasova et al. [22]. The authors obtained a complex site-specific equation for single-level C 2 n , similar to nonML empirical models proposed by, e.g., Sadot [23], Wang et al. [24], Arockia Bazil Raj et al. [25]. The difference between these empirical models and MOST-based models is that they are not derived physically but aim to capture the local optical turbulence strength well. Going beyond single-level C 2 n estimates, Pierzyna et al. [26] proposed a physics-inspired framework to derive a data-driven similarity theory of optical turbulence that combines C 2 n observations from multiple levels within the surface layer into a single non-dimensional model. In contrast to the near-surface studies, Su et al. [27] proposed a method to estimate a vertical C 2 n profile by modeling OT strength as a function of radiosonde measurements. However, all these approaches require in-situ data, rendering them unsuitable for temporal extrapolation as targeted in the present study. Additionally, in-situ-based models need sensors deployed at the site of interest for C 2 n estimation, which hinders their applicability to sites without instrumentation. To alleviate this issue, we proposed OTCliM, which parameterizes C 2 n based on reanalysis data, which are available globally and for multiple decades into the past.

Step 1: Train model on overlapping 1-year time series

Step 2: Make predictions using trained model and full reference dataset

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00520v1-figure-4a578c1b-de4f-4883-846b-10e2d33ed932.jpg)

Step 3: Create climatology from extrapolated

FIG. 1: Proposed OTCliM approach to extrapolate a measured 1-year time series of optical turbulence strength (yellow) to multiple years (orange) based on ERA5 reference data (blue). Robust yearly C 2 n statistics can be obtained from the extrapolated data.

## III. METHODOLOGY

## A. Gradient boosting regression

Our proposed OTCliM approach aims to extrapolate one year of observed C 2 n to multiple years to obtain a robust statistical description of OT strength at a particular site. An overview of OTCliM is given in Fig. 1, which we base on the measure-correlate-predict (MCP) framework popular in wind energy [28, 29]. In the first step, C 2 n is measured for one year (yellow) and then correlated to a reference dataset. For OTCliM, we utilize variables from the ERA5 reanalysis [30] extracted at the location of the C 2 n measurements as reference data and gradient boosting machines (GBM) to regress the ERA5 time series to C 2 n where they overlap (blue). In the second step, the trained model is utilized to extrapolate C 2 n based on ERA5 to multiple years (prediction step in MCP), which can then be used to obtain site-specific seasonal statistics as presented in step 3.

We utilize GBMs for the regression step (cf. § III A) because they are known for their high performance in nonlinear tabular regression tasks. For comparison, we also train GBM models that use in-situ observations instead of ERA5 as input data, which serve as performance baselines for OTCliM. Similarly, the traditional W71 model is evaluated with ERA5 inputs to provide a second baseline (cf. § III B). To ensure our GBM models capture physically reasonable dependencies, we quantify the feature importance assigned to each ERA5 input variable (cf. § III C). Finally, the performance of the trained models is assessed in terms of their temporal extrapolation and geographical generalization capabilities (cf. § III D).

Following the MCP framework presented in Fig. 1, we aim to train one ML model for one site using training data from one year. The corresponding training data are the input data X s,t ∈ R n × p and the target vector y s,t ∈ R n , where s ∈ S is one out of multiple sites S and t ∈ T is one out multiple training years T available at that site. The ML task is to learn a regression function ˆ f s,t that approximates a sample y i ∈ y s,t from the target vector based on a sample x i ∈ X s,t of the input data: ˆ f s,t ( x i ) = ˆ y i ≈ y i . X s,t and y s,t contain a matching number of n samples, i.e., the samples of the one-year time series, where X s,t is composed of p features, i.e., concurrent time series of different meteorological variables, and y s,t contains the scaled log 10 C 2 n .

The OT strength C 2 n is log 10 -transformed because it varies over multiple orders of magnitude throughout the day, which is challenging to capture in non-log space. Also, the range of yearly C 2 n variation varies between sites, so log 10 C 2 n is scaled before training to make the performance scores and feature importance values comparable between sites. More details about the C 2 n preprocessing are presented in § IVA.

The GBM regression models ˆ f s,t are trained using the AutoML library FLAML [31]. FLAML is a timeconstrained hyperparameter optimizer that aims to find the optimal GBM model configuration within a specified time budget. FLAML optimizes not only the hyperparameters of a single algorithm but also explores switching between algorithms, here, the two popular algorithms

XGBoost [32] and LightGBM [33]. The internal loss functions of XGBoost and LightGBM are also treated as hyperparameters to be optimized by FLAML with L 1 -norm, L 2 -norm, and Huber norm as options. The hyperparameter optimization uses 5-fold cross-validation based on the training data and aims to minimize the rootmean-squared error between the predicted and hold-out y . Each model is trained for 45 min on 8 CPU cores of an Intel Xeon 2648R, resulting in 6 core hours per model.

## B. Baseline models

We employ two baseline models to put the performance of OTCliM into perspective. The first model is the traditional empirical W71 model as given by Eqns. 1 and 3. Since the similarity function g ζ ( ) does not contain site-specific model coefficients, the accuracy of the W71 C 2 n estimates is expected to be lower than the sitespecific OTCliM models. Secondly, we utilize GBM models trained with in-situ observations as input instead of reanalysis. Reanalysis datasets are coarse-resolution model outputs (spatial resolution typically O (10 km)), which are expected to miss some local and complex processes and patterns modulating turbulence. Such processes would be contained in the in-situ data, so the performance of in-situ GBM models is expected to be higher than that of ERA5-based OTCliM models. These baselines allow us to disentangle the influence of the modeling approach (traditional vs. GBM) from that of the input data (ERA5 vs. in-situ). Appendix A gives more technical details about both baseline models.

## C. Physical consistency checking using feature importance

To assess if the trained models captured physically consistent relationships and to potentially discover currently unknown dependencies between features and the regression target, we quantify the importance of each feature for the prediction using SHapely Additive exPlanation (SHAP) values [34]. SHAP values explain how a trained model ˆ f s,t arrives at its predictions ˆ on a per-sample y basis. For each sample x i , the prediction ˆ y i is explained as

$$\hat { f } _ { s, t } ( x _ { i } ) = \hat { y } _ { i } = \mathbb { E } \left [ \hat { y } _ { i } \right ] + \sum _ { j = 1 } ^ { p } \phi _ { j, i } \quad \quad ( 4 ) \quad \text{src} \\ \mathcal { T } \, \searrow \, \text{Trh} \,.$$

where ϕ j,i are the SHAP values of the i -th sample and E [ ˆ ] y i is the expected value of all predictions. In plain words, the j -th SHAP value ϕ j,i describes the contribution of the j -th feature to the prediction assuming a local linear model. A large step ϕ j,i compared to the other SHAP values indicates that the j -th feature contributes strongly to that prediction and must be important. Consequently, a global view of the importance of the j -th feature can be obtained by averaging the j -th SHAP values for all n samples of the dataset [35] such that

$$I _ { j } = \frac { 1 } { n } \sum _ { i = 1 } ^ { n } | \phi _ { j, i } |.$$

Note that SHAP values explain a feature's contribution to a model's prediction ˆ , not to the true value y i y i . Accurate feature importance (FI) estimates, thus, require well-performing models ˆ fs, t where ˆ y i closely approximates y i for training and testing data. As predictions typically align more closely with true values in training data than testing data, we compute SHAP values using the former. This approach allows us to examine the model's internal workings without the influence of the generalization error. Moreover, it aligns with practical scenarios where all available data is used for training the final model, meaning that the FI analysis must be performed using the training dataset [35].

## D. Performance evaluation

The performance of the trained OTCliM models is quantified through the Pearson correlation coefficient r and the root-mean-squared error (RMSE) ϵ :

$$r = \frac { \sum _ { i = 1 } ^ { n } ( y _ { i } - \overline { y } ) ( \hat { y } _ { i } - \overline { \hat { y } } ) } { n \, \sigma _ { y } \, \sigma _ { \hat { y } } },$$

$$\epsilon = \sqrt { \frac { 1 } { n } \sum _ { i = 1 } ^ { n } { ( y _ { i } - \hat { y } _ { i } ) } ^ { 2 } }.$$

In Eq. 6, □ and σ □ are the mean and standard deviation of true ( y ) and estimated (ˆ) y scaled C 2 n , respectively. Two validation strategies are employed to assess the temporal and geographical extrapolation capabilities of the OTCliM models. Note that only the testing data used to compute the scores differ between the evaluation strategies - the training data and the trained model remain the same.

a. Temporal extrapolation For the MCP application - the temporal extrapolation -, the model ˆ f s,t trained on data from year t ∈ T needs to accurately predict OT strength at the same site for the held-out testing years T \ t : ˆ f s,t ( x i ) = ˆ y i with x i ∈ X s, T \ t and ˆ y ∈ ˆ y s, T \ t . The scores are computed between the true held-out data y s, T \ t and the predictions ˆ y s, T \ t and are denoted r s,t and ϵ s,t . Training happens in a round-robin fashion, so one model is trained per training year and evaluated for the hold-out years. This process is repeated until all years are used for training once. The average performance of these models for a specific station s , called the MCP score, is obtained by averaging the scores across all training years: □ s = ⟨ □ s,t ⟩ t ∈T .

TABLE I: Flux stations of the New York State Mesonet used to benchmark the OTCliM approach. The three urban stations are marked with (*).

| ID     | Nearest City   |   Latitude, ◦ N |   Longitude, ◦ W | Environment        |
|--------|----------------|-----------------|------------------|--------------------|
| BELL   | Belleville     |           43.79 |            76.11 | Lake               |
| BKLN ∗ | Brooklyn       |           40.63 |            73.95 | Coastal / Urban    |
| BURT   | Burt           |           43.32 |            78.75 | Lake               |
| CHAZ   | Chazy          |           44.9  |            73.46 | Valley             |
| FRED   | Fredonia       |           42.42 |            79.37 | Lake               |
| ONTA   | Ontario        |           43.26 |            77.37 | Lake               |
| OWEG   | Owego          |           42.03 |            76.26 | Plateau            |
| PENN   | Penn Yan       |           42.66 |            76.99 | Mountainous / Lake |
| QUEE ∗ | Queens         |           40.73 |            73.82 | Coastal / Urban    |
| REDF   | Redfield       |           43.62 |            75.88 | Lake               |
| REDH   | Red Hook       |           42    |            73.88 | Valley             |
| SCHU   | Schuylerville  |           43.12 |            73.58 | Valley             |
| SOUT   | Southold       |           41.04 |            72.47 | Coastal            |
| STAT ∗ | Staten Island  |           40.6  |            74.15 | Coastal / Urban    |
| VOOR   | Voorheesville  |           42.65 |            73.98 | Valley             |
| WARS   | Warsaw         |           42.78 |            78.21 | Mountainous / Lake |
| WHIT   | Whitehall      |           43.49 |            73.42 | Valley             |

b. Cross-site evaluation The second evaluation strategy probes the geographical generalization capability of a model trained on station s when applied to another station ˜ s ∈ S \ s . The trained model is tasked to predict OT strength at ˜ given s the full multi-year input data X ˜ s, T at ˜ from which the cross-site s scores using r and RMSE can be obtained: □ ( s,t ) → (˜ s, T ) . For readability, we average the cross-site scores achieved by the individual models per training site s : □ s → ˜ s = ⟨ □ ( s,t ) → (˜ s, T ) ⟩ t ∈T .

## IV. DATASET

To evaluate OTCliM's performance, we train GBM models for 17 diverse locations across New York State (cf. § IVA). Five years of C 2 n are available at each site for which we select collocated and concurrent ERA5 data (cf. § IVB).

## A. New York State Mesonet (NYSM)

The New York State Mesonet (NYSM) comprises 127 standard weather stations spread across New York State, United States of America, and has been fully operational since 2018 [36]. These weather stations measure routine meteorological parameters such as 2 m-temperature, 10 m-wind speed and direction, surface pressure, and several other variables. The sampling rate of these measurements is on the order of seconds, and final values are reported as 10 min aggregates (mean and variance). However, measuring C 2 n requires higher-frequency observations to resolve the inertial range of turbulence. Such high-frequency measurements are available at a subset of 17 stations - the flux stations - additionally equipped with Campbell Scientific CSAT3A sonic anemometers mounted at 9 m height. The NYSM sonic anemometers measure the three wind components and the sonic temperature at f s = 10Hz [36], which is high enough to obtain C 2 n .

We utilize the data from these NYSM flux stations to obtain training targets for the OTCliM models. The flux stations are placed in diverse topographical and climatological environments as listed in Tab. I. The stations BKLN, QUEE, and STAT are located on rooftops in urban environments where measurements are strongly influenced by their immediate surroundings. Neighboring buildings can, for example, cast shadows or cause wakes, which influence radiation, wind, and, therefore, also the local turbulence [37]. Since we expect these urban stations to behave differently than the rural stations, they are marked with (*) in the table. For each of the 17 stations, five years of measurements are available. Following the notation introduced previously, the set of flux stations is denoted S , and T represents the set of five training years 2018 - 2022 available at each site. The corresponding target vector y s,t is obtained by estimating C 2 n from sonic anemometer measurements, applying quality assurance and control (QA/QC) steps, and scaling the OT data to make them comparable between different sites. These stages are described in detail below.

a. Structure function approach As noted earlier, C 2 n quantifies the strength of refractive index fluctuations due to density fluctuations caused by turbulent temperature fluctuations. Therefore, C 2 n can be obtained by measuring the temperature fluctuations quantified by the coefficient of the second-order temperature structure function C 2 T , which is defined in the inertial range as

$$S _ { T } ^ { 2 } ( \Delta x ) = \langle ( T ( x ) - T ( x + \Delta x ) ) ^ { 2 } \rangle \stackrel { \bot } { = } C _ { T } ^ { 2 } \Delta x ^ { 2 / 3 }. \quad ( 8 )$$

The structure function is computed from 5 min nonoverlapping windows of the sonic temperature signal, and C 2 T is obtained by fitting a 2 / 3 slope to the inertial range of S 2 T in log-log coordinates.

b. Quality assurance and control Since our GBM models are purely data-driven, the quality of the fitted model strongly depends on the quality of the training data, so we apply several QA/QC steps to obtain a clean dataset. First, C 2 T has to be determined from the slope of the inertial range, i.e., slopes of 2 / 3. Since the inertial range is not always resolved or well captured in every 5 min window, only C 2 T values with well-fitted slopes are kept (2 / 3 ± 5% slope, R 2 &gt; 0 95). . A similar procedure was followed by He and Basu [38] for the analysis of simulated data. Second, S 2 T in Eq. 8 is defined spatially as a function of separation distance ∆ x , but the sonic anemometer measures T temporally. Taylor's frozen turbulence assumption is commonly made to convert the temporal sonic signal into a spatial form using the horizontal mean wind speed M and ∆ x = M/f s . However, Taylor's assumption might break down in low wind conditions, so all C 2 T values for M &lt; 1 m s -1 are discarded. Finally, precipitation strongly affects sonic temperature measurements [39], so only C 2 T values for dry conditions

can be used. This filter is also justified from an operational perspective because optical telescopes are typically not operated in the rain, and FSOC links are strongly attenuated by snow [40]. While the QA/QC procedure discards ∼ 50% of the data throughout the year, it still contains ca. fifty thousand samples per site and training year spread throughout diverse meteorological conditions, which we view as adequate.

c. Distribution scaling The quality-controlled C 2 T time series is converted to C 2 n using the Gladstone equation (Eq. 2) and then scaled. The scaling is needed because the yearly distribution of OT strength differs between sites. Some locations show stronger OT (shifted distribution) or a larger range of C 2 n values (scaled distribution) than other sites. To make the performance scores and SHAP FI values comparable between sites, we scale the site-specific log 10 C 2 n before training based on the 50% log-range such that

$$\left [ \log _ { 1 0 } C _ { n } ^ { 2 } ( t ) \right ] _ { \text{scaled} } = \frac { \log _ { 1 0 } C _ { n } ^ { 2 } ( t ) - \left [ \log _ { 1 0 } C _ { n } ^ { 2 } \right ] _ { p 2 5 } } { \left [ \log _ { 1 0 } C _ { n } ^ { 2 } \right ] _ { p 7 5 } - \left [ \log _ { 1 0 } C _ { n } ^ { 2 } \right ] _ { p 2 5 } }. \ \ ( 9 ) _ { \text{scaled} }.$$

Here, the subscripts p25 and p75 indicate the 25th and 75th percentiles of the original distribution used for scaling.

## B. ERA5 reanalysis data

The ERA5 reanalysis [30] serves as the input dataset for OTCliM from which C 2 n shall be estimated. ERA5 is available globally from 1950 until the present and, thus, includes our NYSM training locations and times. We extract ERA5 input data for all NYSM stations and training years at the ERA5 grid points closest to the respective stations. ERA5 grid points represent 1/4-degree grid boxes (ca. 30 km by 30 km), so by selecting the closest point, we obtain data representative for the grid box containing the station. Figure 2 depicts the location of the flux stations (grey markers) and their collocated grid boxes (orange square) and also gives an impression of how the coarse ERA5 represents (a) land, sea, and (b) terrain.

We aim to incorporate all variables linked to the two processes driving atmospheric turbulence: wind shear and buoyancy. Commonly used variables include sensible heat flux and friction velocity, as utilized in W71. However, an advantage of using ML over deriving analytical equations from theory is the ability to include variables that may indirectly influence turbulence. For example, the ERA5 gravity wave dissipation rate (GWD) could be significant in complex mountainous regions where orographic gravity wave drag is known to modulate momentum fluxes [41, 42]. If there is a relationship between GWD and C 2 n , ML models will likely identify and utilize it, potentially revealing new dependencies.

A complete list of the ERA5 variables selected as features is presented in Tab. II. Many listed features

## (a) Land sea mask

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00520v1-figure-52c2e0d8-3df8-4099-a00f-cb4e2806f643.jpg)

- (b) Geopotential height

FIG. 2: ERA5 representation of NYSM domain with true locations of NYSM flux stations (grey) and corresponding ERA5 grid box (orange) containing the stations. Urban sites are marked with (*).

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00520v1-figure-d700302a-1ccc-48fc-864b-c5df79dd689e.jpg)

are (partially) redundant and/or (partially) correlated. Since it is initially unknown which features are most suitable for estimating near-surface C 2 n , we aim for a broad feature space. During training, the GBM algorithm will identify and base its predictions on the most important features. Through the post-training feature importance analysis, we can identify which features the trained model deems relevant and assess their physical relevance. It is usually preferred in ML to minimize the number of features and reduce model complexity. However, for treebased algorithms like GBM used in this study, irrelevant features are automatically mostly disregarded, mitigating this concern.

The features in Tab II without ERA5 variable names are so-called engineered features, meaning that they are derived from one or more ERA5 variables. In the following, we detail the variable selection and the feature engineering.

a. Shear-related features Table II lists multiple wind-related features aimed at capturing wind shear, i.e.,

TABLE II: ERA5 variables and variables derived thereof serve as features for the OTCliM approach. Derived features do not have ERA5 variable names and are marked with (-).

| Category            | Symbol        | Grouping   | ERA5 variable name   | Description                                                                                        | Unit    | Unit    |
|---------------------|---------------|------------|----------------------|----------------------------------------------------------------------------------------------------|---------|---------|
| Wind                | M 100         |            | -                    | 100m horizontal wind speed, see eq. 10                                                             | ms      | - 1     |
| Wind                | M 10          |            | -                    | 10m horizontal wind speed, see eq. 10                                                              | ms      | - 1     |
| Wind                | G 10          |            | i10fg                | Instantaneous 10m wind gust                                                                        | ms      | - 1     |
| Shear               | u ∗           |            | zust                 | Friction velocity                                                                                  | ms      | - 1     |
| Shear               | α             |            | -                    | Shear coefficient of power-law wind profile, see eq. 11                                            | -       |         |
| Shear               | β             |            | -                    | Directional shear, | X 10 - X 100 |                                                                | ◦       |         |
| Fetch               | cos ( X 100 ) | Σ X        | -                    | Cosine of 100m wind direction, see eq. 12                                                          | -       |         |
| Fetch               | cos ( X 10 )  | Σ X        | -                    | Cosine of 10m wind direction, see eq. 12                                                           | -       |         |
| Fetch               | sin ( X 100 ) | Σ X        | -                    | Sine of 100m wind direction, see eq. 12                                                            | -       |         |
| Fetch               | sin ( X 10 )  | Σ X        | -                    | Sine of 10m wind direction, see eq. 12                                                             | -       |         |
| Buoyancy            | Q H           |            | ishf                 | Instantaneous surface sensible heat flux                                                           | Wm      | - 2     |
| Temperature         | T 2           |            | t2m                  | 2m temperature                                                                                     | K       |         |
| Temperature         | T 0           |            | skt                  | Skin temperature                                                                                   | K       |         |
| Temperature         | T soil        |            | stl1                 | Soil temperature                                                                                   | K       |         |
| Stability           | ∆ T 0         | Σ∆ T       | -                    | Temperature difference T soil - T                                                                  | K       |         |
| Stability           | ∆ T 1         | Σ∆ T       | -                    | 0 Temperature difference T 0 - T 2                                                                 | K       | - 2     |
| Radiation           | R S ↓ , diff  | Σ R        | msdwswrf             | Mean diffuse downwelling short-wave radiation flux at surface                                      | Wm      | - 2     |
| Radiation           | R S ↓ , dir R | Σ R        | msdrswrf             | Mean direct downwelling short-wave radiation flux at surface                                       | Wm      |         |
| Radiation Radiation | S R L         | Σ R        | msnswrf msdwlwrf     | Mean net short-wave radiation flux at surface Mean downwelling long-wave radiation flux at surface | Wm Wm   | - 2 - 2 |
|                     | ↓             | Σ R        | msnlwrf              | Mean net long-wave radiation flux at surface                                                       | Wm      | - 2     |
| Radiation           | R L           | Σ R Σcc    | lcc                  | Low cloud cover fraction (0 - 1)                                                                   |         |         |
| Cloud Cloud         | lcc tcc       | Σcc        | tcc                  | Total cloud cover fraction (0 - 1)                                                                 | - -     |         |
| Moisture            | ∆ T d         |            | -                    | Dew point spread, T - T 2                                                                          | K       |         |
|                     |               |            |                      | d 2                                                                                                |         |         |
| Moisture Moisture   | T d 2 Q       |            | d2m ie               | 2m dew point temperature Instantaneous moisture flux                                               | K kgm - | 2 s - 1 |
| Moisture            | q Q L         |            | mslhf                | Mean surface latent heat flux                                                                      | Wm      | - 2     |
| Auxilary            |               |            | blh                  | Boundary layer height                                                                              | m       |         |
| Auxilary            | h i BLD       |            | bld                  | Boundary layer dissipation rate                                                                    | Jm      | - 2     |
| Auxilary            | GWD           |            | gwd                  | Gravity wave dissipation                                                                           | Jm -    | 2       |
| Auxilary            | CAPE          |            | cape                 | Convective available potential energy                                                              | Jkg -   | 1       |
| Auxilary            | P 0           |            | msl                  | Mean sea level pressure                                                                            | Pa      |         |
|                     |               |            | dndza                | inside trapping layer                                                                              | m -     | 1       |
| Auxilary            | ∂n/∂z a       |            |                      | Mean vertical gradient of refractivity                                                             |         |         |
| Auxilary            | h db          |            | dctb                 | Duct base height                                                                                   | m       |         |
| Auxilary Auxilary   | h tpb         |            | tplb                 | Trapping layer base height                                                                         | m       |         |
|                     | h tpt         |            | tplt                 | Trapping layer top height                                                                          | m       |         |
| Time                | cos hr ′      | Σhr ′      | -                    | Cosine of normalized hour of the day                                                               | -       |         |
| Time                | sin hr ′      | Σhr ′      | -                    | Sine of normalized hour of the day                                                                 | -       |         |
| Time                | cos day ′     | Σday ′     | -                    | Cos of normalized day of the year                                                                  | -       |         |
| Time                | sin day ′     | Σday ′     | -                    | Sin of normalized day of the year                                                                  | -       |         |
| Time                | cos month ′   | Σmonth ′   | -                    | Cos of normalized month of the year                                                                | -       |         |
| Time                | sin month ′   | Σmonth ′   | -                    | Sin of normalized month of the year                                                                | -       |         |

the vertical change of the wind. Like traditional C 2 n parameterizations, such as W71, we select the friction velocity u ∗ as the first shear-related feature. We also include the horizontal wind speed obtained from the zonal and meridional wind fields of ERA5 U z and V z - as

$$M _ { z } = \sqrt { U _ { z } ^ { 2 } + V _ { z } ^ { 2 } } \quad \quad ( 1 0 ) \quad \text{Tr}$$

at z = 10m and z = 100m above ground. Assuming a power-law wind profile of the form M z ( ) = M ref ( z/z ref ) α between the 10 m and the 100 m level, we utilize the exponent of the power law α as an additional shear-related feature:

$$\alpha ( z _ { 1 }, z _ { 2 } ) = \frac { \log M _ { z _ { 2 } } - \log M _ { z _ { 1 } } } { \log z _ { 2 } - \log z _ { 1 } }. \quad \ \ ( 1 1 ) \quad \ \ f f { \cdot } { \cdot }$$

The directional shear, i.e., wind turning with height, is included through the absolute angular difference between 10 m and 100 m wind direction defined as β = | X 10 -X 100 | where the wind direction X z is given as

$$X _ { z } = \arctan 2 \left ( - U _ { z }, - V _ { z } \right ). \quad \quad \left ( 1 2 \right )$$

The wind direction also serves as a proxy for upstream effects, or fetch, that might influence the observed turbulence. For example, atmospheric turbulence measured at a station close to the coast can be very different if the wind blows from land to sea or sea to land. The periodicity of X z is accounted for by including the sines and cosines of X z as fetch features instead of X z directly.

b. Buyoancy and stability The prime candidate to reflect the influence of buoyancy on C 2 n is the sensi-

ble heat flux Q H . Sensible heat flux is also featured in W71 because it captures buoyancy and static atmospheric stability (in ERA5's convention, Q H &gt; 0 during stable/nighttime conditions and Q H &lt; 0 during unstable/daytime conditions). Additionally, static stability can be estimated from temperature gradients, so we utilize the absolute values of 2 m-temperature, skin temperature, and soil temperature, as well as the differences between them as complementary features.

c. Surface Energy Budget The strength of buoyancy depends on the radiation that reaches the surface, so we also include ERA5 radiation fluxes and cloud cover as features to complement the buoyancy information from Q H . The simplified surface energy budget describes how radiation forces the surface fluxes. If the small soil heat flux into the ground is neglected, the steady-state energy balance at the surface is [43]

$$Q _ { L } + Q _ { H } \approx R _ { S } + R _ { L }, \text{ \quad \ \ } ( 1 3 ) \text{ \ \ } \text{raw}$$

where Q L is the latent heat flux due to evaporation and R S and R L are the net shortwave and longwave surface radiation fluxes. The net fluxes R □ = R □ ↓ -R □ ↑ are the difference between incoming/downwelling (index ↓ ) and the reflected/upwelling radiation (index ↑ ). ERA5 contains the net and downward radiation fluxes as variables. Further, it splits the shortwave downwelling radiation R S ↓ into a direct component R S , ↓ dir and a diffuse component R S , ↓ diff resulting from scattering by clouds. Since clouds affect how much shortwave radiation reaches the surface and how much longwave radiation from the surface is reflected, we include ERA5's low cloud cover ( z &lt; ca. 2 km) and total cloud cover into the dataset. Close to water bodies, Q L is expected to be higher due to stronger evaporation leading to lower Q H , and, thus, lower C 2 n . Therefore, Q L and other moisture-related ERA5 variables are listed as moisture-related features in Tab. II.

d. Auxiliary features Finally, we include several auxiliary features more loosely related to wind shear and buoyancy, such as boundary layer dissipation rate (BLD), convective available potential energy (CAPE), or the aforementioned gravity dissipation rate (GWD). Also, certain daily and seasonal patterns exist in meteorology, which we aim to capture through synthetic time-dependent features. Based on the timestamp of the data point, we compute the sines and cosines (for periodicity) of the normalized hour of the day (hr ′ = 2 π hr / 24 h), the normalized day of the year (day ′ = 2 π day 365d), / and the normalized month of the year (month ′ = 2 π month 12mo). /

## V. RESULTS: FEATURE IMPORTANCE

The physics governing optical turbulence are the same everywhere - in urban or rural environments or for inland stations or coastal sites. The modulating processes are always buoyancy and wind shear. However, the

ERA5 features that best reflect and predict these processes locally can differ between sites. To assess and quantify potential differences, we present the SHAPbased feature importance values for all trained models in Fig. 3. To make the FI analysis less verbose, we make use of the linearity of the SHAP values (cf. Eq. 4) and group the SHAP values of features related to similar physical processes as presented in Tab. II. These groups are prefixed with Σ.

Panel (a) depicts the importance distributions of all features for all models. This global view reveals that buoyancy-related features such as radiation Σ R , sensible heat flux Q H , and cloud cover Σ cc are key variables for C 2 n prediction. Also, the shear exponent α is picked up, indicating that the OTCliM models identify wind shear as a modulating factor of turbulence. The dependency on wind direction Σ X suggests that upstream effects influence C 2 n prediction at many sites, while atmospheric stability is likely captured next to Q H by the temperature differences Σ∆ T and the boundary layer height h i . All of the above aligns with our physical knowledge of atmospheric turbulence, indicating that the models picked up physically meaningful relations between ERA5-generated variables and C 2 n . However, the boxplots also show outliers for several features, indicating that not all models agree with the global view.

The FI values are displayed per station in Fig. 3(b), allowing us to assess the outliers in more detail. The ten features with the highest average FI for each location are highlighted in orange. Since different stations exhibit different top ten features, 20 features are shown in total, but their markers are colored in light grey if the features are not in the top ten of the site in question. For reference, the black lines for each feature represent the global average FI values according to panel (a). Comparing all per-station plots on a high level reveals two key points. First, although the FI values associated with the five models trained for each site exhibit some scatter, the overall importance of features at one site is consistent between models. Therefore, we conclude that the models captured features representative of that site throughout the training years. Secondly, the color coding highlights that different features are important for different sites, suggesting that the processes dominating turbulence vary between locations.

The three urban stations, BKLN, QUEE, and STAT, and the coastal SOUT site deviate most obviously from the mean FI distribution. In particular, BKLN and QUEE show an above-average dependency of estimated C 2 n on wind direction. Both stations are located on rooftops and have some tall buildings in their neighborhood [36]. Such inhomogeneous urban conditions are known to influence local measurements [37], so it seems reasonable that turbulence strength observed at BKLN and QUEE depends more strongly on wind direction compared to other more rural locations. All four stations (urban + SOUT) also show below-average dependency on Q H , which we view as a feature capturing the

FIG. 3: SHAP value-based feature importance (FI) of all OTCliM models (a) aggregated and (b) per training site. In (b), urban stations are marked with (*), the top 10 features of each station are highlighted in orange, and the global FI averages are indicated as black dashes. Panel (c) repeats the geopotential height from Fig. 2 to aid geographical interpretation of the results.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00520v1-figure-c3a86d3c-4bd6-4b08-bc06-46f2346e5885.jpg)

diurnal cycle. Instead, the models seem to have picked Σhr or ′ h i as predictors for this information. This behavior differs from all other models, where Q H is typically the second most important input, in line with more traditional physics-based C 2 n models, such as W71. This shift, however, need not be viewed as unfavorable but as a demonstration of the power of ML-based modeling to shift from traditional to 'unconventional' features if complex flow conditions require it. The downside is that we expect such models to perform more poorly when applied to other sites where the more traditional features are relevant. Especially, BKLN is expected to generalize poorly because it also depends strongly on the temperature features T □ , deemed irrelevant by almost all other models.

Less drastic but distinct differences are also visible between the non-urban models. Wind shear α , for example, has lower than average importance assigned for models trained on lake shores (FRED, BURT, ONTA). The models of another set of stations (WHIT, VOOR, SCHU, PENN, OWEG) picked up the gravity-wave dissipation rate (GWD). These sites are located in valleys or mountainous areas where gravity waves could modulate near-surface turbulence [41, 42]. Still, the dependency is small, and stations CHAZ and REDH, located at the end of valleys but still surrounded by mountains, do not depend on GWD. A more detailed study of the sites' climatologies would be needed, which is beyond the scope of this work.

Overall, the FI values of most models represent the known physical dependencies of atmospheric turbulence. That supports our confidence that our OTCliM approach is well suited for MCP. The 'unconventional' features picked up by some models are viewed as an advantage of ML-based methods to arrive at accurate predictions even in complex environments. However, we assume that geographical generalization will be more difficult for such models, which is addressed later in this paper.

## VI. RESULTS: OTCLIM PERFORMANCE

After establishing that all OTCliM models picked up physically meaningful dependencies, we turn our analysis to quantifying their performance. The foundation for this analysis are 85 models trained individually for the 17 NYSMstations and the five training years. The temporal extrapolation capabilities of each model are quantified in § VIA to assess the suitability of OTCliM for MCP. To put the MCP scores achieved by our approach into perspective, we compare them to the two baseline models, W71 and the in-situ-based GBM models. Since we have a network of stations available, we also assess the geographical generalizability of the OTCliM models by applying models trained on one site to all other sites. The results of this cross-site evaluation study are discussed in § VIB.

## A. Temporal extrapolation

The temporal extrapolation performance of each model ˆ f s,t with respect to the correlation coefficient r and the scaled RMSE ϵ is presented in panels (a) and (b) of Fig. 4, respectively. The heatmaps for both metrics show the MCP scores □ s,t with their site-averaged values □ s being compared to the baseline models in the accompanying scatter plots. The heatmap reveals similar patterns for both metrics: the score variance across training years (rows) is low for almost all stations, while the site-specific performance (columns) varies notably. In other words, regardless of the year of C 2 n observations used for training, the MCP performance stays consistent, whereas C 2 n at some locations (e.g., QUEE, BKLN, SOUT) is easier to predict from ERA5 than at others (e.g., VOOR, WHIT). It is unexpected but impressive that the urban models marked with (*) are the highest performing ones due to the typically complex urban climatology, which makes traditional modeling difficult [44]. Panels (c) and (d) show a few randomly drawn batches of observed C 2 n with their corresponding predictions. These curves highlight that the overall performance at SOUT (c) is better than that of WHIT (d) for two reasons. First, the observed C 2 n (black) at WHIT exhibits more high-frequency oscillations than SOUT, which are missing in the predictions (red). Second, complex short-duration events seem to occur more frequently for WHIT, which are only partially captured. Both effects indicated that WHIT is subject to more complex flows than SOUT and that OTCliM misses these details. We attribute this to missing details in ERA5, as will be discussed later. For now, we conclude that the models for both locations capture the general trends of C 2 n well and that the observed intrasite performance variation is primarily due to smoothed predictions. For practical applications such as obtaining a C 2 n climatology, the large-scale behavior of C 2 n is the most relevant, so the smoothing is of limited concern.

The insight that the training year selection has a low influence on the model performance has high practical relevance. It suggests that models can be trained on archive data and do not necessarily require recent observations iff the site's climatology remains comparable between the training period and the envisioned period of model application. Drastic changes in the station's surroundings, e.g., through construction projects or climate change, would likely break this assumption and, thus, the model's applicability.

The variability across models trained for different stations is likely due to ERA5 not fully capturing all local effects modulating the turbulence. For example, some landscape features, such as small lakes around PENN and WARS, are smaller than the ERA5 resolution, so they are missed. The relevance of these local features is highlighted when the OTCliM performance (black circles) is compared to that of the in-situ baseline models in the r s and ϵ s scatter plots. While the in-situ models exhibit a small performance scatter around their mean accuracy

## (a) Correlation coefficient r (higher is better)

## (b) Root-mean-squared error ϵ (lower is better)

FIG. 4: Performance of OTCliM models compared to the baseline models. The heatmaps show the performance of individual OTCliM models trained with one year of C 2 n observations from site s when evaluated on the four hold-out years of the same site. The scatter plots compare the site-averaged performance against the baseline estimates. Seven-day batches of C 2 n are randomly drawn to compare observations (black) against predictions (red) for a high-accuracy (c) and a lower-accuracy OTCliM model (d), with complex turbulence conditions shaded in grey.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00520v1-figure-6a53363d-61b6-4f68-9f57-1cb0cf23d461.jpg)

(dashed lines), a generally widening performance gap is visible, which sets off small (QUEE) and grows toward the lower end of the performance ranking. Remember that only the input data differs between OTCliM and insitu (ERA5 vs. in-situ observations); the GBM training framework remains the same. Consequently, the growing gap between the model variants can be attributed to ERA5 holding less predictive power for some locations than others. Expectedly, the models with lower accuracy (e.g., WHIT, VOOR, PENN) are located close to lakes or in mountainous areas, which are usually subject to complex flows that ERA5 might not well capture. Surprisingly, the ERA5-based urban models perform almost on par with the in-situ-based models. This result is im- pressive, given that urban flows are complex and likely not well represented in ERA5. Nevertheless, as discussed earlier, GBM seemingly circumvented this issue by shifting from traditional to unconventional features in urban cases. These features still seem to hold enough predictive power for the high performance observed here.

A final comparison between OTCliM and the traditional W71 parameterization in the scatter plots of panels (a) and (b) shows that OTCliM clearly outperforms W71. That is because the underlying similarity function (cf. Eq. 3) is empirically determined based on the flow over a flat, unobstructed plain [14] and does not adapt to local topography or climatology. This limitation becomes especially evident for the urban sites QUEE and BKLN

/negationslash

/negationslash

∈ S

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00520v1-figure-c07aea61-4ef5-4a52-a1a3-23da26be6aaa.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00520v1-figure-ff5113ff-96a6-4e12-84fb-5755c0cd7c03.jpg)

∈ S

(b) Relative and absolute root-mean-squared error ϵ (lower is better)

/negationslash

/negationslash

∈ S

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00520v1-figure-f184b7fb-0bb4-43bb-b8e4-9b8cb4af6246.jpg)

∈ S

FIG. 5: Cross-site evaluation performance of OTCliM models. The row of a heatmap presents the performance of the model trained on site s when evaluated on all other sites ˜ s ∈ S . The left heatmaps show the performance degradation compared to the MCP case, i.e, the relative scores □ ′ s → ˜ s = □ s → ˜ s / □ s → s , and the right ones show the absolute scores for reference. The histograms depict the distributions of both heatmaps.

where r s is much below the W71 average (dashed line). Our OTCliM models, on the other hand, perform well at these locations, highlighting the advantage of accounting for local climatology in complex atmospheric conditions. In summary, the presented scores demonstrate the capability of our OTCliM approach to accurately extrapolate a 1-year C 2 n time series to multiple years with ERA5 input for a diverse set of locations.

## B. Cross-site evaluation

Next, we assess the geographical generalization capabilities of the OTCliM models by assessing their performance when being evaluated across different sites. Each model ˆ f s,t trained on one site s is asked for C 2 n predictions based on five years of ERA5 input data of all other sites (˜ s ∈ S \ s ). The resulting cross-site (c/s) correlation coefficient and RMSE scores are presented in Fig. 5. Each cell ( s, s ˜) of the heatmaps corresponds to the average score achieved by the five OTCliM models trained on site s (rows) and evaluated on site ˜ (columns). s The heatmaps on the left present the scores normalized by the MCP scores, reflecting the performance degradation of a model when applied to other sites compared to performance on the hold-out years of the original training site. The heatmaps on the right show the absolute scores, and the histograms in the middle display the distribution of both relative and absolute scores with diagonals excluded, i.e., without the MCP performance. The rows and columns of all heatmaps are ordered by mean generalization performance quantified as ⟨ r ′ s → ˜ s ⟩ ˜ s ∈S\ s .

Overall, the geographical generalization works very well for the non-urban sites. For half of all cases, correla-

→

→

tion and RMSE performance degrades by not more than 16% for r and 25% for ϵ compared to the original MCP scores. Considering 75% of the cases, corresponding to the c/s scores of almost all non-urban sites, the absolute performance falls as low as r = 0 56 and . ϵ = 0 63, which . is still significantly better than the W71 baseline scores of ⟨ r s ⟩ = 0 3 and . ⟨ ϵ s ⟩ = 1 15 (cf. Fig 4). . In other words, any non-urban OTCliM model applied to any other nonurban site yields significantly better performance than the traditional approach.

As expected from the FI analysis, the urban models do not perform well in non-urban locations. Both relative and absolute c/s scores in Fig. 5 are low for STAT, QUEE, and BKLN, resulting in long tails in the c/s histograms. However, the urban models perform better at other urban sites, as indicated by the small square of higher performance in the bottom right corner of the c/s matrices. It seems that the performance of urban generalization could be linked to the degree of urbanization. BKLN is the most urban site, QUEE has more vegetation than BKLN, and STAT has more vegetation than QUEE, which is reflected in the progression of c/s performance: BKLN generalizes most poorly and depends on the most uncommon features, QUEE performs slightly better, and STAT performs even better with FI values similar to the average case. Nevertheless, this observation might be unique to these New York City stations, which are all relatively close to each other and experience similar mesoscale conditions.

## VII. CONCLUSION

This study presents OTCliM, a gradient boostingbased measure-correlate-predict approach to obtain climatologies of optical turbulence strength from one year of C 2 n observations and multi-year ERA5 reference data. A feature importance (FI) analysis based on SHAP values revealed that most models, especially the non-urban ones, learn similar dependencies. The dominating features reflect the known and expected physical dependencies, suggesting that the OTCliM models captured physically meaningful relations between ERA5 and C 2 n . The FI results also suggest features such as boundary layer dissipation rate or gravity wave dissipation rate as C 2 n dependencies, which are not considered in traditional models. This result demonstrates how machine learning models could hint at new insights into turbulent processes. In contrast to most stations, the three urban locations differed significantly from the average FI distribution. We assume that such models learned very different inputoutput relations, hindering them from generalizing well to other sites. In a temporal extrapolation study, we demonstrate that the OTCliM models accurately predict four held-out years of C 2 n from ERA5 inputs across 17 diverse sites in New York State. The models' performances stay consistent regardless of the year of data selected for training, suggesting archive data can also be used for training. The trained models are also found suitable for stations located in complex environments such as valleys or cities, which are often difficult to model traditionally. Minor accuracy variations between models from different sites are visible but can be attributed to local details missing in the coarse ERA5 input data. For some regions in the world, regional reanalysis datasets are available at higher resolution (e.g., Europe [45] or Australia [46]) or with more assimilated data (e.g., North America [47]), which might help to reduce this problem. Already with ERA5 input, OTCliM models significantly outperform traditional analytical C 2 n models and achieve scores close to those of GBM models trained on in-situ inputs. A geographical generalization study also showed that applying non-urban models to other non-urban sites yields high scores in many cases, indicating that these models generalize well. The performance only degrades significantly when the urban models are evaluated at non-urban sites. While these models performed very well in the temporal extrapolation test case, the low scores for cross-site application indicate that the urban models learned very site-specific relations.

The key conclusion for practice from our work is that one year of C 2 n observations is sufficient to obtain reliable site-specific C 2 n climatologies using OTCliM. Measuring C 2 n at sites of interest for only one year instead of multiple years would free up instruments sooner. As a result, site survey costs could be reduced, or more locations could be surveyed for the same costs. OTCliM's high performance in urban environments is highly relevant for FSOC, where terminals can be expected to be located in both rural and urban locations. The good geographical generalization results of the non-urban models indicate that future OTCliM models trained on a few stations can potentially be used to estimate the spatial distribution of C 2 n over larger areas. Such C 2 n maps could provide a first indication of locations suitable for FSOC or astronomy with respect to optical turbulence. However, the region-specific features identified in the FI analysis for lake or mountain environments also suggest that the models must account for local effects. Developing such multi-region models requires more work before large-scale C 2 n estimates can be considered reliable. As a first step, users can perform an FI analysis on their trained models and use the results as guidance for expected geographical generalization. For example, if the models pick up unconventional features, they will likely not generalize to other sites, but if features are more traditional, they might.

Finally, we believe that OTCliM can be highly relevant also for C 2 n forecasts. The rapidly advancing ML weather prediction (MLWP) models, such as Graphcast [48] or AIFS [49], are typically trained on ERA5 and, therefore, produce ERA5-like forecasts. Currently, not all variables required for OTCliM are available from MLWP, but in principle, OTCliM can be used to translate the MLWP forecasts to C 2 n forecasts. Similarly, one could train OTCliM on historical data from the Global Forecasting System (GFS) [50], a traditional numerical weather predic-

tion system, and convert the GFS forecasts to C 2 n forecasts.

Besides all these advantages, the critical assumption of OTCliM is that the observed year of C 2 n represents the temporal extrapolation range. If the surroundings change drastically through, e.g., construction, this assumption breaks down, and OTCliM predictions can no longer be considered valid. Also, our proposed approach currently only predicts near-surface C 2 n , but for astronomy and FSOC, higher-level C 2 n or full profiles are also relevant. In a previous study, we presented an ML-based framework that can combine C 2 n from multiple levels into one physics-inspired model [26] for multi-level C 2 n predictions. How this approach can be integrated into OTCliM remains open for future work. In summary, our presented OTCliM approach is shown to be highly accurate in diverse meteorological conditions with the potential for geographical generalization. We believe that OTCliM is a relevant tool for the optical turbulence community for future site surveys and related studies.

## ACKNOWLEDGMENTS

This publication is part of the project FREE - Optical Wireless Superhighways: Free photons (at home and in space) (with project number P19-13) of the research programme TTW-Perspectief, which is (partly) financed by the Dutch Research Council (NWO). This research is made possible by the New York State (NYS) Mesonet. Original funding for the NYS Mesonet (NYSM) buildup was provided by Federal Emergency Management Agency grant FEMA-4085-DR-NY. The continued operation and maintenance of the NYSM is supported by National Mesonet Program, University at Albany, Federal and private grants, and others.

## Appendix A: Details about baseline models

This section details how the traditional W71 C 2 n parameterization and the in-situ-based GBM models are setup and utilized as performance baselines for our proposed OTCliM approach.

## 1. W71 C 2 n parameterization

All variables on which the W71 equations (cf. Eqns. 1 - 3) depend are available from ERA5 (cf. Tab. II). Only

[1] J. W. Hardy, Adaptive Optics for Astronomical Telescopes , Vol. 16 (Oxford University Press, USA, 1998).

[2] H. Kaushal and G. Kaddoum, Optical Communication in Space: Challenges and Mitigation Tech-

the dynamical sensible heat flux Q H from ERA5 needs to be converted to its kinematic form w θ ′ ′ as w θ ′ ′ = -Q / ρc H ( p ), where ρ and c p are the density and specific heat capacity of air [43]. Also, the sign of Q H is flipped because ERA5 treats an upward Q H as negative, whereas W71 expects the upward w θ ′ ′ to be positive. After this conversion, we utilize the same ERA5 data extracted for the 17 NYSM stations that we use to train the OTCliM models to evaluate the W71 equations. The results are 17 five-year C 2 n time series, which can be compared to the observed C 2 n evolution at each station. The corresponding RMSE and r scores are presented in Fig. 4.

## 2. In-situ GBM models

These upper-bound models are similar to the OTCliM models but use different input data. The OTCliM model of a specific site s is a GBM model (cf. III A) that employs ERA5 input data extracted from the grid box containing site s and C 2 n observed at site s . Instead of ERA5, one can utilize in-situ weather data observed at s . The resulting GBM in-situ models are quite similar, for example, to work presented by Wang and Basu [18], Jellen et al. [21], or Pierzyna et al. [26]. As for selecting ERA5 variables (cf. IVB), we aim to capture buoyancy and wind shear with multiple features based on the instruments deployed at the NYSM stations. As buoyancy features, we compute the temperature gradient Γ = ( T 9 -T 2 ) / (9 m -2 m) from observed 2 m and 9 m temperature, and utilize the observed 30 min sensible heat flux and observed incoming radiation. Assuming the velocity at the ground to be close to zero, we obtain a crude estimate of the bulk wind shear S = M / 10 10 m where M 10 is the horizontal windspeed at 10 m. Additionally, observed friction velocity u ∗ , latent heat flux w q ′ ′ , and the dew point spread ∆ T d = T 2 -T d 2 are included as features, where T d 2 is the 2 m dew point temperature. To account for atmospheric stability, Γ and S are combined into a bulk Richardson number as Ri = ( g/T ) Γ /S 2 with g = 9 81 m s . -2 . Upstream effects are captured by adding sine and cosine of the 10 m wind direction. All these variables form the insitu dataset, which is used to train GBM models in the same fashion as OTCliM: each in-situ model is trained on one year of observations from site s and evaluated on the hold-out years of the same site. The resulting scores are presented as upper baseline in Fig. 4.

niques, IEEE Commun. Surv. Tutorials 19 , 57 (2017), arXiv:1705.10630 [cs, math].

[3] A. Jahid, M. H. Alsharif, and T. J. Hall, A contemporary survey on free space optical communication: Poten-

- tials, technical challenges, recent advances and research direction, Journal of Network and Computer Applications 200 , 103311 (2022).
- [4] F. Hill, J. Beckers, P. Brandt, J. Briggs, T. Brown, W. Brown, M. Collados, C. Denker, S. Fletcher, S. Hegwer, T. Horst, M. Komsa, J. Kuhn, A. Lecinski, H. Lin, S. Oncley, M. Penn, R. Radick, T. Rimmele, H. SocasNavarro, and K. Streander, Site testing for the Advanced Technology Solar Telescope, in SPIE Astronomical Telescopes + Instrumentation , edited by L. M. Stepp (Orlando, Florida , USA, 2006) p. 62671T.
- [5] M. Sch¨ ock, S. Els, R. Riddle, W. Skidmore, T. Travouillon, R. Blum, E. Bustos, G. Chanan, S. G. Djorgovski, P. Gillett, B. Gregory, J. Nelson, A. Ot´rola, J. Seguel, a J. Vasquez, A. Walker, D. Walker, and L. Wang, Thirty Meter Telescope Site Testing I: Overview, PASP 121 , 384 (2009).
- [6] C. Fuchs and F. Moll, Ground Station Network Optimization for Space-to-Ground Optical Communication Links, J. Opt. Commun. Netw. 7 , 1148 (2015).
- [7] S. Poulenard, M. Crosnier, and A. Rissons, Ground Segment Design for Broadband Geostationary Satellite With Optical Feeder Link, J. Opt. Commun. Netw. 7 , 325 (2015).
- [8] T. V. Pham, H. Yamano, and I. Susumu, A placement method of ground stations for optical satellite communications considering cloud attenuation, IEICE Commun. Express 12 , 568 (2023).
- [9] E. Masciadri, J. Vernin, and P. Bougeault, 3D mapping of optical turbulence using an atmospheric numerical model - I. A useful tool for the ground-based astronomy, Astron. Astrophys. Suppl. Ser. 137 , 185 (1999).
- [10] M. Pierzyna, R. Saathof, and S. Basu, A multi-physics ensemble modeling framework for reliable C 2 n estimation, in Environmental Effects on Light Propagation and Adaptive Systems VI , Vol. 12731 (SPIE, 2023) pp. 185-191.
- [11] M. Pierzyna, O. Hartogensis, S. Basu, and R. Saathof, Intercomparison of flux, gradient, and variance-based optical turbulence ( C 2 n ) parameterizations, Appl. Opt. 10.1364/AO.519942 (2024).
- [12] F. Beyrich, O. K. Hartogensis, H. A. de Bruin, and H. C. Ward, Scintillometers, in Springer Handbook of Atmospheric Measurements , Springer Handbooks, edited by T. Foken (Springer International Publishing, Cham, 2021) pp. 969-997.
- [13] M. Beason, G. Potvin, D. Sprung, J. McCrae, and S. Gladysz, Comparative analysis of C 2 n estimation methods for sonic anemometer data, Appl. Opt. 63 , E94 (2024).
- [14] J. C. Wyngaard, Y. Izumi, and S. A. Collins, Behavior of the Refractive-Index-Structure Parameter near the Ground, J. Opt. Soc. Am. 61 , 1646 (1971).
- [15] A. S. Monin and A. M. Obukhov, Basic laws of turbulent mixing in the surface layer of the atmosphere, Contrib. Geophys. Inst. Acad. Sci. 151 (1954).
- [16] L. C. Andrews and R. L. Phillips, Optical Turbulence in the Atmosphere, in Laser Beam Propagation Through Random Media (SPIE, 2005) pp. 57-83.
- [17] M. J. Savage, Estimation of evaporation using a dualbeam surface layer scintillometer and component energy balance measurements, Agricultural and Forest Meteorology 149 , 501 (2009).
- [18] Y. Wang and S. Basu, Using an artificial neural network approach to estimate surface-layer optical turbulence at
- Mauna Loa, Hawaii, Opt. Lett. 41 , 2334 (2016).
- [19] A. M. Vorontsov, M. A. Vorontsov, G. A. Filimonov, and E. Polnau, Atmospheric Turbulence Study with Deep Machine Learning of Intensity Scintillation Patterns, Applied Sciences 10 , 8136 (2020).
- [20] C. Jellen, J. Burkhardt, C. Brownell, and C. Nelson, Machine learning informed predictor importance measures of environmental parameters in maritime optical turbulence, Appl. Opt. 59 , 6379 (2020).
- [21] C. Jellen, M. Oakley, C. Nelson, J. Burkhardt, and C. Brownell, Machine-learning informed macrometeorological models for the near-maritime environment, Appl. Opt. 60 , 2938 (2021).
- [22] L. A. Bolbasova, A. A. Andrakhanov, and A. Y. Shikhovtsev, The application of machine learning to predictions of optical turbulence in the surface layer at Baikal Astrophysical Observatory, Monthly Notices of the Royal Astronomical Society 504 , 6008 (2021).
- [23] D. Sadot, Forecasting optical turbulence strength on the basis of macroscale meteorology and aerosols: Models and validation, Opt. Eng. 31 , 200 (1992).
- [24] H. Wang, B. Li, X. Wu, C. Liu, Z. Hu, and P. Xu, Prediction model of atmospheric refractive index structure parameter in coastal area, Journal of Modern Optics 62 , 1336 (2015).
- [25] A. Arockia Bazil Raj, J. Arputha Vijaya Selvi, and S. Durairaj, Comparison of different models for ground-level atmospheric turbulence strength ( C 2 n ) prediction with a new model according to local weather data for FSO applications, Appl. Opt. 54 , 802 (2015).
- [26] M. Pierzyna, R. Saathof, and S. Basu, Π-ML: A dimensional analysis-based machine learning parameterization of optical turbulence in the atmospheric surface layer, Opt. Lett. 48 , 4484 (2023).
- [27] C. Su, X. Wu, S. Wu, Q. Yang, Y. Han, C. Qing, T. Luo, and Y. Liu, In situ measurements and neural network analysis of the profiles of optical turbulence over the Tibetan Plateau, Monthly Notices of the Royal Astronomical Society 506 , 3430 (2021).
- [28] J. A. Carta, S. Vel´zquez, and P. Cabrera, A review of a measure-correlate-predict (MCP) methods used to estimate long-term wind characteristics at a target site, Renewable and Sustainable Energy Reviews 27 , 362 (2013).
- [29] S. Kartal, S. Basu, and S. J. Watson, A decision-treebased measure-correlate-predict approach for peak wind gust estimation from a global reanalysis dataset, Wind Energy Science 8 , 1533 (2023).
- [30] H. Hersbach, B. Bell, P. Berrisford, S. Hirahara, A. Hor´ anyi, J. Mu˜ noz-Sabater, J. Nicolas, C. Peubey, R. Radu, D. Schepers, A. Simmons, C. Soci, S. Abdalla, X. Abellan, G. Balsamo, P. Bechtold, G. Biavati, J. Bidlot, M. Bonavita, G. Chiara, P. Dahlgren, D. Dee, M. Diamantakis, R. Dragani, J. Flemming, R. Forbes, M. Fuentes, A. Geer, L. Haimberger, S. Healy, R. J. Hogan, E. H´ olm, M. Janiskov´, S. Keeley, P. Laloyaux, a P. Lopez, C. Lupu, G. Radnoti, P. Rosnay, I. Rozum, F. Vamborg, S. Villaume, and J.-N. Th´paut, The ERA5 e global reanalysis, Q.J.R. Meteorol. Soc. 146 , 1999 (2020).
- [31] C. Wang, Q. Wu, M. Weimer, and E. Zhu, FLAML: A Fast and Lightweight AutoML Library, in Proceedings of Machine Learning and Systems , Vol. 3 (2021) pp. 434447.
- [32] T. Chen and C. Guestrin, XGBoost: A Scalable Tree Boosting System, in Proceedings of the 22nd ACM

- SIGKDD International Conference on Knowledge Discovery and Data Mining , KDD '16 (Association for Computing Machinery, New York, NY, USA, 2016) pp. 785794.
- [33] G. Ke, Q. Meng, T. Finley, T. Wang, W. Chen, W. Ma, Q. Ye, and T.-Y. Liu, LightGBM: A Highly Efficient Gradient Boosting Decision Tree, in Advances in Neural Information Processing Systems , Vol. 30 (Curran Associates, Inc., 2017).
- [34] S. M. Lundberg and S.-I. Lee, A unified approach to interpreting model predictions, in Proceedings of the 31st International Conference on Neural Information Processing Systems , NIPS'17 (Curran Associates Inc., Red Hook, NY, USA, 2017) pp. 4768-4777.
- [35] C. Molnar, Interpretable Machine Learning: A Guide for Making Black Box Models Explainable , second edition ed. (Christoph Molnar, Munich, Germany, 2022).
- [36] J. A. Brotzge, J. Wang, C. D. Thorncroft, E. Joseph, N. Bain, N. Bassill, N. Farruggio, J. M. Freedman, K. Hemker, D. Johnston, E. Kane, S. McKim, S. D. Miller, J. R. Minder, P. Naple, S. Perez, J. J. Schwab, M. J. Schwab, and J. Sicker, A Technical Overview of the New York State Mesonet Standard Network, Journal of Atmospheric and Oceanic Technology 37 , 1827 (2020).
- [37] World Meteorological Organization, Volume III -Observing Systems , Guide WMO-No. 8 (WMO, Geneva, 2023).
- [38] P. He and S. Basu, Development of similarity relationships for energy dissipation rate and temperature structure parameter in stably stratified flows: A direct numerical simulation approach, Environ Fluid Mech 16 , 373 (2016).
- [39] R. Zhang, J. Huang, X. Wang, J. A. Zhang, and F. Huang, Effects of precipitation on sonic anemometer measurements of turbulent fluxes in the atmospheric surface layer, J. Ocean Univ. China 15 , 389 (2016).
- [40] ITU, Prediction Methods Required for the Design of Terrestrial Free-Space Optical Links , Recommendation Rec. ITU-R P.1814 (International Telecommunication Union, 2007).
- [41] D. K. Lilly, Wave Momentum Flux-A GARP Problem, Bulletin of the American Meteorological Society 53 , 17 (1972).
- [42] T. N. Palmer, G. J. Shutts, and R. Swinbank, Alleviation of a systematic westerly bias in general circula-
- tion and numerical weather prediction models through an orographic gravity wave drag parametrization, Quart J Royal Meteoro Soc 112 , 1001 (1986).
- [43] R. B. Stull, An Introduction to Boundary Layer Meteorology (Kluwer Academic Publishers, Dordrecht, 1988).
- [44] M. W. Rotach, R. Vogt, C. Bernhofer, E. Batchvarova, A. Christen, A. Clappier, B. Feddersen, S.-E. Gryning, G. Martucci, H. Mayer, V. Mitev, T. R. Oke, E. Parlow, H. Richner, M. Roth, Y.-A. Roulet, D. Ruffieux, J. A. Salmond, M. Schatzmann, and J. A. Voogt, BUBBLE an Urban Boundary Layer Meteorology Project, Theor. Appl. Climatol. 81 , 231 (2005).
- [45] Copernicus Climate Change Service, CERRA sub-daily regional reanalysis data for Europe on single levels from 1984 to present (2021).
- [46] C.-H. Su, N. Eizenberg, P. Steinle, D. Jakob, P. FoxHughes, C. J. White, S. Rennie, C. Franklin, I. Dharssi, and H. Zhu, BARRA v1.0: The Bureau of Meteorology Atmospheric high-resolution Regional Reanalysis for Australia, Geoscientific Model Development 12 , 2049 (2019).
- [47] F. Mesinger, G. DiMego, E. Kalnay, K. Mitchell, P. C. Shafran, W. Ebisuzaki, D. Jovi´, J. Woollen, E. Rogers, c E. H. Berbery, M. B. Ek, Y. Fan, R. Grumbine, W. Higgins, H. Li, Y. Lin, G. Manikin, D. Parrish, and W. Shi, North American Regional Reanalysis, Bulletin of the American Meteorological Society 87 , 343 (2006).
- [48] R. Lam, A. Sanchez-Gonzalez, M. Willson, P. Wirnsberger, M. Fortunato, F. Alet, S. Ravuri, T. Ewalds, Z. Eaton-Rosen, W. Hu, A. Merose, S. Hoyer, G. Holland, O. Vinyals, J. Stott, A. Pritzel, S. Mohamed, and P. Battaglia, Learning skillful medium-range global weather forecasting, Science 0 , eadi2336 (2023).
- [49] S. Lang, M. Alexe, M. Chantry, J. Dramsch, F. Pinault, B. Raoult, M. C. A. Clare, C. Lessig, M. MaierGerber, L. Magnusson, Z. B. Bouall`gue, A. P. Nemesio, e P. D. Dueben, A. Brown, F. Pappenberger, and F. Rabier, AIFS - ECMWF's data-driven forecasting system (2024), arXiv:2406.01465 [physics].
- [50] National Centers for Environmental Prediction, NCEP GFS 0.25 Degree Global Forecast Grids Historical Archive (2015).
Maximilian Pierzyna Sukanta Basu Rudolf Saathof
OTCliM: generating a near-surface climatology of optical turbulence strength ($C_n^2$) using gradient boosting
This study introduces OTCliM (Optical Turbulence Climatology using Machine
learning), a novel approach for deriving comprehensive climatologies of
atmospheric optical turbulence strength ($C_n^2$) using gradient boosting
machines. OTCliM addresses the challenge of efficiently obtaining reliable
site-specific $C_n^2$ climatologies near the surface, crucial for ground-based
astronomy and free-space optical communication. Using gradient boosting
machines and global reanalysis data, OTCliM extrapolates one year of measured
$C_n^2$ into a multi-year time series. We assess OTCliM's performance using
$C_n^2$ data from 17 diverse stations in New York State, evaluating temporal
extrapolation capabilities and geographical generalization. Our results
demonstrate accurate predictions of four held-out years of $C_n^2$ across
various sites, including complex urban environments, outperforming traditional
analytical models. Non-urban models also show good geographical generalization
compared to urban models, which capture non-general site-specific dependencies.
A feature importance analysis confirms the physical consistency of the trained
models. It also indicates the potential to uncover new insights into the
physical processes governing $C_n^2$ from data. OTCliM's ability to derive
reliable $C_n^2$ climatologies from just one year of observations can
potentially reduce resources required for future site surveys or enable studies
for additional sites with the same resources.