Impact of a Turbulent Ocean Surface on Laser Beam Propagation

The roughness of the ocean surface significantly impacts air-to-sea imaging, oceanographic monitoring, and optical communication. Most current and previous methods for addressing this roughness and its impact on optical propagation are either entirely statistical or theoretical, or are ‘mixed methods’ based on a combination of statistical models and parametric-based physical models. In this paper, we performed experiments in a 50-foot-wave tank on wind-generated waves, in which we varied the wind speed to measure how the surface waves affect the laser beam propagation and develop a geometrical optical model to measure and analyze the refraction angle and slope angle of the laser beam under various environmental conditions. The study results show that the laser beam deviations/distortions and laser beam footprint size are strongly related to wind speed and laser beam incidence angle.


Introduction
The ocean covers more than 70 percent of the world's surface. Over the years, people have explored the oceans to satisfy basic scientific curiosity and more pragmatic concerns such as shipping routes and schedules, oil field maintenance, or tactical surveillance. Nevertheless, according to The National Oceanic and Atmospheric Administration (NOAA), about 95% of the oceans remain unexplored [1,2]. To continue monitoring the ocean and marine activities, autonomous underwater vehicles (AUV), remotely operated underwater vehicles (ROV), and other sensors have been deployed to gather, collect, and transmit data about this environment.
The demand for reliable, high-speed communication is accelerating; it is estimated [3] that by 2030 there will be more than 40 billion devices connected to the internet simultaneously. These devices include computers, smart devices, lidars, and, eventually, autonomous underwater vehicles. The concept of the Internet of Things (IoT) was invented in 1985 [4]; in 2012, the Internet of Underwater Things (IoUT) was first discussed. The IoUT is defined as "the network of smart interconnected underwater objects". The smart objects could be several types of underwater sensors, autonomous underwater vehicles (AUVs), autonomous surface vehicles (ASVs), buoys, ships, etc.
Furthermore, with the rapid development of autonomous underwater vehicles (AUVs) and unmanned aerial vehicles (UAVs), there is an increasing need for reliable communication between these two types of platforms for efficient information exchange. In recent years, considerable work has been directed toward underwater acoustic, RF, and wireless optical technologies. Despite the advances in underwater acoustic communication (UAC), it suffers a significant latency (slow speed of sound in water), which leads to considerable Inter-Symbol Interference (ISI). Thus, UAC links have extremely limited data rates and do If the air-sea interface were quiescent, then it would be a specular surface with two layers with different refractive indices. However, the air-sea interface is rough due to ocean waves [8], especially wind waves. According to the Beaufort wind scale, wind can be divided into 12 scales or forces [9]. In this paper, our upper limit is Beaufort 3 (wind speed is below 5.5 m/s) because, at this level, the laboratory waves start breaking, which leads to gas bubbles and whitecap generation. As a result, the performance of the optical communication path across the ocean surface will be degraded or even hard to establish [10].
Ocean surface waves are typically measured by two approaches [11]. The traditional method uses in situ buoys, such as those employed by the National Data Buoy Center (NDBC) [12]. A second method is to take instantaneous images of the ocean waves, then process and analyze them to obtain the oceanic wave spectrum. The wave spectrum provides the distribution of wave energy or variance contained over different temporal frequencies and spatial wavelengths on the ocean surface. It can be described by either a onedimensional (1D) frequency spectrum E(f) or by a two-dimensional (2D) frequency If the air-sea interface were quiescent, then it would be a specular surface with two layers with different refractive indices. However, the air-sea interface is rough due to ocean waves [8], especially wind waves. According to the Beaufort wind scale, wind can be divided into 12 scales or forces [9]. In this paper, our upper limit is Beaufort 3 (wind speed is below 5.5 m/s) because, at this level, the laboratory waves start breaking, which leads to gas bubbles and whitecap generation. As a result, the performance of the optical communication path across the ocean surface will be degraded or even hard to establish [10].
Ocean surface waves are typically measured by two approaches [11]. The traditional method uses in situ buoys, such as those employed by the National Data Buoy Center (NDBC) [12]. A second method is to take instantaneous images of the ocean waves, then process and analyze them to obtain the oceanic wave spectrum. The wave spectrum provides the distribution of wave energy or variance contained over different temporal frequencies and spatial wavelengths on the ocean surface. It can be described by either a one-dimensional (1D) frequency spectrum E(f) or by a two-dimensional (2D) frequency directional spectrum E(f,θ) [13]. A review of ocean wave spectra is provided by Ryabkova [14], and a summary of widely used models is given in Section 2.2 below.
In this paper, we present experiments on wind-generated waves, in which we varied the wind speed to measure how the surface waves affect the optical propagation of a laser beam. We considered laser beams with two different diameters from below the glassbottomed wave tank. Data from specific points on the surface were captured by a highspeed camera above the wave tank. These measurements were compared to established models of ocean wave spectra. The main objective of this work is to understand the ∇ 2 ϕ = 0 f or x ∈ R, y ∈ R, z ∈ (−h, η(x, y, t)). ( The kinematic boundary condition that defines the free surface is: The dynamic boundary at the free surface is that the curvature balances the pressure jump due to surface tension so that: where T is the kinematic surface tension, and g is the acceleration of gravity. The boundary condition at the bottom is that there is no flow through the bottom. The vertical velocity vanishes there so that: The solution of the linearized version of this boundary-value problem, along with (1), is that the velocity potential is: where κ = √ κ 2 + l 2 is the 2D wavenumber. The corresponding wave phase speed is and, in general, may depend on the wavenumber. Several classifications of water waves are obtained from Equation (6).

1.
The first classification is with respect to the depth through the size of κh. For κh 1, tanh(κh) ≈1, the effects of depth are neglected, and the result is "deep-water waves", also called "short waves". For these waves, C p is a function of wavenumber, so waves with different wavelengths travel at different speeds. Thus, deep-water waves are dispersive waves. The waves discussed in this paper are deep-water waves. For κh 1, tanh(κh) ≈ κh, and the results are "shallow-water waves", also called "long-waves". For these waves, C p is approximately independent of wavenumber. Thus, shallowwater waves are approximately non-dispersive. For κh =1, there is no approximation on tanh(κh). Such waves are dispersive.

2.
The second classification is with respect to the Bond number [17]: which measures the relative importance of capillary forces versus gravitational forces. For β = 1 (using T = 73.0 cm 3 s 2 and g = 981 cm s 2 ), the wavelength is a critical value of 1.71 cm, and the two restoring forces balance. This wavelength corresponds to the minimum phase speed of about C p = 23.1 cm/s for deep-water waves. For β 1, the wavelengths are shorter than the critical value, and capillary forces dominate. The dispersion relation for these capillary waves is well approximated by: For β 1, the wavelengths are longer than the critical value, and gravitational forces dominate. The dispersion relation for these gravity waves is well approximated by:

3.
A third classification takes into account weak nonlinearity. In particular, Equation (6) holds when the wave slope is very small, aκ→ 0. If one allows for finite but weak nonlinearity so that aκ 1, one finds that capillary-gravity waves and capillary waves may spread energy spectrally through resonant triad and quartet interactions as well as modulational instabilities [18]. Gravity waves on finite depth or deep water spread energy through modulational instabilities and resonant quartet interactions. Thus, even in the absence of wind, these instabilities and interactions may cause complicated two-dimensional surface patterns [19] from freely propagating waves.

4.
A fourth classification is with respect to the presence or absence of wind-forcing. Waves are classified as being either "sea" or "swell", where seas are the waves that feel the influence of wind-forcing, and swelling are the waves that have propagated away from the influence of the wind. Because of the wind-forcing, seas are steeper than swells; they have a larger value of aκ than the swells. Because deep-water waves are dispersive, the swells sorted themselves into narrow-banded spectra, with the longer waves traveling faster than the shorter waves. Their frequencies are smaller, and their wavelengths are longer than those of the sea. A typical separation frequency In summary, the waves considered herein are seas comprised of deep-water, capillarygravity waves of finite amplitude. In addition, they are fetch-limited. Fetch is the distance over which the wind blows. The minimum duration time, t min is the time for waves to travel from the beginning of the fetch to a distance, r. If the wind blows for a time larger than t min , then the wave height at the position, r, stops growing, and the waves there reach (on average) a steady-state height [21]. For the experiments herein, measurements are obtained for times, t > t min .

Models of Ocean Wave Spectra
Waves on the ocean surface are not monochromatic and are typically modeled as a superposition. The wind puts energy into a broad spectrum of wave modes so that one may write the surface displacement in terms of a Fourier Transform: where k and l are the x and y wavenumbers, ω is the corresponding frequency, A(k, l, ω) is the Fourier amplitude, and one adds the complex conjugate to obtain real values for the free surface displacement. Such a representation is too complicated because one needs complete information about the storm that generated the waves to solve for A. In addition, this view does not account for wave damping primarily due to wave breaking and energy transfer due to nonlinear interactions during propagation. Instead, investigators have developed a wide variety of models for spectra [14]. Pierson-Moskowitz [22] used similarity theory to find a closed-form representation of the frequency spectrum for fully-developed seas, a situation in which the energy input by wind and the dissipation due to breaking are balanced. The Pierson-Moskowitz energy spectral density is where α = 8.1 × 10 −3 , δ = 0.74, ω 0 = g/U, and U is the wind speed measured above the surface, typically at 10 m, if possible. This speed, U, is used since the friction velocity at the interface is not measured. However, data obtained during the Joint North Sea Wave Observation Project (JONSWAP) [23], showed that the ocean wave spectrum is typically not fully developed. To account for nonlinear interactions and the fetch-dependent balance of energy sources, sinks, and energy transfer, the authors in [21] modified the Pierson-Moskowitz spectrum and used data from the project to develop the JONSWAP spectrum, which varies with respect to U (measured 10 m above the surface in reference to JONSWAP paper) and x, the fetch. Here: where ω p is the frequency at the spectrum's peak, and g is the acceleration of gravity. The inclusion of x allows the spectrum to vary with distance from the wind source, as per the data obtained in [24]. A summary of several significant wave models and spectra is listed in Table 1. In Section 5.2, Figure 14, we compare measured spectra obtained from laser measurements to those obtained using the PM and JONSWAP spectra. Table 1. A summary of candidate models for oceanic waves.

Spectrum/Model Description References
Gerstner Waves (1802) It is based on Navier-Stokes equation by describing a particle's motion on the surface as a circular motion to provide an approximate to simulate the air-water interface. [12,25] Phillips (1954) A fully developed sea is considered deep water. It is widely used in real-time simulation of oceanic waves. [9,25,26] Neumann (1955) It is valid for only fully developed sea, and it is valid for only gravity waves regime. [8]

Pierson-Moskowitz (1964)
It is based on the Phillips equilibrium range representing a fully developed sea. It is designed to describe gravity waves over infinite fetch. [13,24] JONSWAP (1964) A modified Pierson-Moskowitz spectrum with enhanced peak and fetch dependent factors. It is valid only for limited fetch and infinite water depth. [13,27,28] TMA (1985) Developed as an extension of the JONSWAP spectrum for finite water depth. [26,29] Majumdar & Brown (1992) Probabilistic method applied to investigate the influence of the wavy air-sea interface on the laser beam transmission based on the Gram-Charlier model. [30] Apel (1994) Modified version of the JONSWAP spectrum includes improved capillary and gravity-capillary wave predictions. It is developed for shallow water with short fetch winds (100-1000) m. [10,31,32] Elfouhaily (1997) Using data observations from previous models, a unified directional spectrum for long and short wind-driven waves based on the Apel wave spectrum. [33]

Laser Propagation at the Air-Sea Interface
Numerous researchers have carried out theoretical and experimental investigations of oceanic surface waves' influence on laser propagation. The relationship between sea surface conditions and the accuracy of airborne LiDAR bathymetry (ALB) was investigated at [10,34,35]. In [1], they experimentally demonstrated high-speed optical wireless communication for both the uplink and downlink by employing an OFDM transmission of 32-QAM and single-mode pigtailed green-light laser diode (LD). Although a data rate of 5.5 Gbps was achieved over a 26 m air-water link with accurate pointing between the transmitter and the receiver, their experiment assumed the oceanic surface was static.
Wang et al. [36] investigated the turbulent propagation of radial partially coherent beams and Gaussian Schell model beams in an air-sea hybrid link scenario while ignoring the turbulent interface. A fast analysis method was proposed to compute the transmittance and the refraction angles through a wavey interface in different wind directions and speeds by simulating the wind-driven dynamic waves for lidar application [37]. The impact of water height rather than wavy interface was investigated at [38]. In [39], they demonstrated a diffuse-line-of-sight communication link through a wavy interface by UV LED as a signal carrier. Majmudar [30] proposed a probabilistic model for the refraction angle of optical propagation at the random air-water interface. A high-speed system for direct optical communications across a water-air interface in a real environment was designed and tested in [40]. AmphiLight was presented and proposed in [41] to enable bidirectional link air-sea interface. Adib [42] proposed hybrid acoustic-RF wireless communication through the water-air interface. For simplicity, the triangular wave facet model was presented and simulated in [43][44][45] using MC ray tracing to compute the sea surface optical reflectance and transmittance after intersecting the interface.
Zhang [46] investigated the effect of large sea surface scale facets on EM scattering by using a capillary wave modification facet scattering model. However, to our knowledge, theoretical and experimental investigations for characterizing the micro-oceanic facets and their influence on laser beam propagation through the interface have not been reported previously.
When a laser beam propagates from underwater to the atmosphere will change direction and likely beam shape due to refraction, as illustrated in Figure 2. The Law of Refraction (Snell's Law) explains the relationship between the incident angle and the refraction angle when the light passes through different media, such that, Irregular (random) ocean surface waves, which are generated by winds, are the major factor in changing beam direction beyond the smooth surface deviation seen in Equation (14) after passing the air-sea interface [47].
speeds by simulating the wind-driven dynamic waves for lidar application [37]. The im pact of water height rather than wavy interface was investigated at [38]. In [39], they demonstrated a diffuse-line-of-sight communication link through a wavy interface by UV LED as a signal carrier. Majmudar [30] proposed a probabilistic model for the refraction angle of optical propagation at the random air-water interface. A high-speed system fo direct optical communications across a water-air interface in a real environment was de signed and tested in [40]. AmphiLight was presented and proposed in [41] to enable bidi rectional link air-sea interface. Adib [42] proposed hybrid acoustic-RF wireless commu nication through the water-air interface. For simplicity, the triangular wave facet mode was presented and simulated in [43][44][45] using MC ray tracing to compute the sea surface optical reflectance and transmittance after intersecting the interface. Zhang [46] investigated the effect of large sea surface scale facets on EM scattering by using a capillary wave modification facet scattering model. However, to ou knowledge, theoretical and experimental investigations for characterizing the micro-oce anic facets and their influence on laser beam propagation through the interface have no been reported previously.
When a laser beam propagates from underwater to the atmosphere will change di rection and likely beam shape due to refraction, as illustrated in Figure 2. The Law o Refraction (Snell's Law) explains the relationship between the incident angle and the re fraction angle when the light passes through different media, such that, = = (14 Irregular (random) ocean surface waves, which are generated by winds, are the major fac tor in changing beam direction beyond the smooth surface deviation seen in Equation (14 after passing the air-sea interface [47].

Experimental Apparatus and Procedures
The experimental apparatus comprised a wave tank, a fan assembly, a laser assembly, and a photographic assembly. Figure 3 shows a photograph of the wave tank. It is 50 ft long, 10 in wide, and can be filled up to 12 in.

Experimental Apparatus and Procedures
The experimental apparatus comprised a wave tank, a fan assembly, a laser assembly, and a photographic assembly. Figure 3 shows a photograph of the wave tank. It is 50 ft long, 10 in wide, and can be filled up to 12 in. We filled it with tap water to a depth of 20 cm for these experiments. The tank has precisely aligned glass sidewalls and a bottom supported by a steel structure. The wind tunnel fits into the wave tank, with its vertical position adjustable to be the desired height (about 1cm) above the water surface. It used a Can Max Fan Mixed Flow Inline Fan with a 10-inch duct diameter that blows up to 1023 ft 3 /min and rotates at a speed up to 2990 rpm. The fan blows into an enclosed Plexiglas chamber the same width as the tank and is 25 cm high and 25 cm long. The downstream wall is a gate with an array of holes. The roof has a slit with a door. We closed the gate and opened the slit so that when the wind was turned on, it blew the air up toward the ceiling, not over the water surface. The fan ramped up to the desired wind speed over about 10 s. After it reached its steady state, we closed the slit and opened the gate simultaneously to create an impulsive start for the wind over the water surface. The air blew from the chamber through three layers of filter material, through a honeycomb of tubes 2.5 in long, then over the water. The airspeed was measured using a TSI, 8465-12 anemometer, and the wind speed value, U, was obtained at a desired height above the water surface for all experiments. While the goal was to investigate how a wavy air-sea interface affects laser propagation, the experiments also provided data on wind-generated waves.
To obtain visual data on wind-generated waves, a light source and a Photron FAST-CAM Mini UX UX100 high-speed camera were placed above the wave tank at a specific angle in the "front" and "back" of the wave direction, respectively. Figure 4 shows a schematic of the setup for these experiments. We used three wind speeds (maximum at 3.3 m/s, medium at 1.9 m/s, and minimum at 0.9 m/s). We took instantaneous images as a video stream of a portion of the air-sea interface, shown as the dashed square in Figure 4, for five to seven seconds at speeds of 1 kf/s and 4 kf/s. The images were compared to those of the flat interface (no wind) to determine the displacement of the air-sea interface.
A doubled Nd:YAG laser operating at 532 nm wavelength was placed below the wave tank with a beam expander to study the effects of laser propagation through a wavy interface. After the expansion process, the beam was reflected by a mirror to obtain the desired location and angle at the receiver. A high-speed digital camera was located at the top of the wave tank to capture the laser beam footprint on a translucent screen 27 cm above the air-water interface. Figure 5 shows a schematic of the set-up for these experiments. We filled it with tap water to a depth of 20 cm for these experiments. The tank has precisely aligned glass sidewalls and a bottom supported by a steel structure. The wind tunnel fits into the wave tank, with its vertical position adjustable to be the desired height (about 1 cm) above the water surface. It used a Can Max Fan Mixed Flow Inline Fan with a 10-inch duct diameter that blows up to 1023 ft 3 /min and rotates at a speed up to 2990 rpm. The fan blows into an enclosed Plexiglas chamber the same width as the tank and is 25 cm high and 25 cm long. The downstream wall is a gate with an array of holes. The roof has a slit with a door. We closed the gate and opened the slit so that when the wind was turned on, it blew the air up toward the ceiling, not over the water surface. The fan ramped up to the desired wind speed over about 10 s. After it reached its steady state, we closed the slit and opened the gate simultaneously to create an impulsive start for the wind over the water surface. The air blew from the chamber through three layers of filter material, through a honeycomb of tubes 2.5 in long, then over the water. The airspeed was measured using a TSI, 8465-12 anemometer, and the wind speed value, U, was obtained at a desired height above the water surface for all experiments. While the goal was to investigate how a wavy air-sea interface affects laser propagation, the experiments also provided data on wind-generated waves.
To obtain visual data on wind-generated waves, a light source and a Photron FAST-CAM Mini UX UX100 high-speed camera were placed above the wave tank at a specific angle in the "front" and "back" of the wave direction, respectively. Figure 4 shows a schematic of the setup for these experiments. We used three wind speeds (maximum at 3.3 m/s, medium at 1.9 m/s, and minimum at 0.9 m/s). We took instantaneous images as a video stream of a portion of the air-sea interface, shown as the dashed square in Figure 4, for five to seven seconds at speeds of 1 kf/s and 4 kf/s. The images were compared to those of the flat interface (no wind) to determine the displacement of the air-sea interface.
A doubled Nd:YAG laser operating at 532 nm wavelength was placed below the wave tank with a beam expander to study the effects of laser propagation through a wavy interface. After the expansion process, the beam was reflected by a mirror to obtain the desired location and angle at the receiver. A high-speed digital camera was located at the top of the wave tank to capture the laser beam footprint on a translucent screen 27 cm above the air-water interface. Figure 5 shows a schematic of the set-up for these experiments.   To analyze the laser propagation through the random interface, a geometric optical model is considered and shown in Figure 6. Let be the incident angle underwater, be the refraction angle at a quiescent interface, be the deviated refraction angle at time( ). The deviations of these angles in the along-wind and cross-wind (assuming independence) directions can be computed using right triangles as follows. The right triangle in the xz-plane with perpendicular sides and ℎ give the angle in the along-wind direction as: Similarly, the right triangle in the yz-plane (out of the page) with perpendicular sides and ℎ give the angle in the cross-wind direction as: where and represent the change in the refraction angle from that of the quiescent interface in the along-wind and cross-wind directions, respectively.   To analyze the laser propagation through the random interface, a geometric optical model is considered and shown in Figure 6. Let be the incident angle underwater, be the refraction angle at a quiescent interface, be the deviated refraction angle at time( ). The deviations of these angles in the along-wind and cross-wind (assuming independence) directions can be computed using right triangles as follows. The right triangle in the xz-plane with perpendicular sides and ℎ give the angle in the along-wind direction as: Similarly, the right triangle in the yz-plane (out of the page) with perpendicular sides and ℎ give the angle in the cross-wind direction as: where and represent the change in the refraction angle from that of the quiescent interface in the along-wind and cross-wind directions, respectively. To analyze the laser propagation through the random interface, a geometric optical model is considered and shown in Figure 6. Let θ i be the incident angle underwater, θ r f be the refraction angle at a quiescent interface, θ r d be the deviated refraction angle at time (t i ). The deviations of these angles in the along-wind and cross-wind (assuming independence) directions can be computed using right triangles as follows. The right triangle in the xz-plane with perpendicular sides di X and h give the angle in the along-wind direction as: Similarly, the right triangle in the yz-plane (out of the page) with perpendicular sides di Y and h give the angle in the cross-wind direction as: where θ rdx and θ rdy represent the change in the refraction angle from that of the quiescent interface in the along-wind and cross-wind directions, respectively. Sensors 2022, 22, x FOR PEER REVIEW 10 of 20 Figure 6. Schematic diagram of the light path through a rough ocean surface. , represent the laser beam centroid at the quiescent interface and at frame , respectively. The laser receiving screen is placed at an approximated mean height ℎ = 27 cm. and are the incident and refraction angles at the quiescent interface case. The deviation in the refraction angle in the along-wind direction depends on the slope angle averaged over the beam footprint at the turbulent interface.
The terms and are the differences between the laser beam centroid "center of mass" at frame and the centroid for the reference frame along the x-axis (along-wind direction) and y-axis (cross-wind direction), respectively. The ℎ value represents the height of the target screen. Then, the total refraction angle in the along-wind direction can be written as: = (17) Similarly, for the cross-wind direction, the total refraction angle is: Snell's law Equation (14) must be modified to include the interfacial slope and so one may determine its angle, as shown in Figure 6. We can then rewrite Equation (14) in the following form: Using Equation (19), the slope angle is calculated as:

Results and Discussion
Here we present data on the two sets of experiments. For the first, we measured the interface displacement from the high-speed camera above the water in front of the wave direction and the diffuse light source in the back of the wave direction (and vice versa) versus three different wind speeds. In the second, we used a laser beam with an expander placed below the wave tank at two different incident angles at the air-water interface (normal incidence = 0° and off-normal incidence = 32°) and the camera placed above the tank, facing 90° downward to capture the laser footprint projected on the receiving screen placed above the surface. The terms di X and di Y are the differences between the laser beam centroid "center of mass" at frame i and the centroid for the reference frame along the x-axis (along-wind direction) and y-axis (cross-wind direction), respectively. The h value represents the height of the target screen. Then, the total refraction angle in the along-wind direction can be written as: Similarly, for the cross-wind direction, the total refraction angle is: Snell's law Equation (14) must be modified to include the interfacial slope θ s and so one may determine its angle, as shown in Figure 6. We can then rewrite Equation (14) in the following form: n w sin(θ i + θ s ) = n a sin(θ r + θ s ) Using Equation (19), the slope angle is calculated as:

Results and Discussion
Here we present data on the two sets of experiments. For the first, we measured the interface displacement from the high-speed camera above the water in front of the wave direction and the diffuse light source in the back of the wave direction (and vice versa) versus three different wind speeds. In the second, we used a laser beam with an expander placed below the wave tank at two different incident angles at the air-water interface (normal incidence θ i = 0 • and off-normal incidence θ i = 32 • ) and the camera placed above the tank, facing 90 • downward to capture the laser footprint projected on the receiving screen placed above the surface.

Water Surface Spectrum
Since oceanic waves and our laboratory waves have random/stochastic, but not necessarily isotropic motions, the statistical measurements of the air-sea interface are different under differing experimental configurations. Therefore, the illumination of the water facets depends on the incident angle, θ i , of the diffused light source and the refraction from the interface into the camera to obtain the wave spectrum [25]. Thus, locating the high-speed camera above the wave tank with wave direction (forward) produces different spectra from when the camera is located in the opposite direction (backward). Figure 7 shows snapshots of the air-sea interface for both cases at the three wind speeds. ward) produces different spectra from when the camera is located in the opposite direction (backward). Figure 7 shows snapshots of the air-sea interface for both cases at the three wind speeds.
Increasing wind speed increases the surface roughness scale, which creates more oceanic surface elements (facets). Each of these individual facets (depending on their number and size) has a slope and tilting angle that is important to understand when the laser propagates and refracts through the air-sea interface [9,48]. In Figure 7a,b, the white arrow indicates the wind direction. In Figure 7b,d,f, the dashed lines represent example facet sizes at a minimum, medium, and maximum wind speed, respectively. The average facet size was measured to be approximately 8.3 cm 2 , 4.9 cm 2 , and 2.1 cm 2 for minimum, medium, and maximum wind speeds, respectively. Figure  8 compares the intensity distribution histogram of the pixel intensity values for the backward case at the three wind speeds that were shown in Figure 7b,d,f. From the histogram, the abscissa presents the grey level intensity, which graphically displays 256 numbers showing the distribution of pixels amongst those greyscale values where 0 means the image is purely black, and 255 means the image is overexposed. The ordinate shows how many pixels there are in each case. Increasing wind speed increases the surface roughness scale, which creates more oceanic surface elements (facets). Each of these individual facets (depending on their number and size) has a slope and tilting angle that is important to understand when the laser propagates and refracts through the air-sea interface [9,48]. In Figure 7a,b, the white arrow indicates the wind direction. In Figure 7b,d,f, the dashed lines represent example facet sizes at a minimum, medium, and maximum wind speed, respectively. The average facet size was measured to be approximately 8.3 cm 2 , 4.9 cm 2 , and 2.1 cm 2 for minimum, medium, and maximum wind speeds, respectively. Figure 8 compares the intensity distribution histogram of the pixel intensity values for the backward case at the three wind speeds that were shown in Figure 7b,d,f. From the histogram, the abscissa presents the grey level intensity, which graphically displays 256 numbers showing the distribution of pixels amongst those greyscale values where 0 means the image is purely black, and 255 means the image is overexposed. The ordinate shows how many pixels there are in each case. tributed in a very tight range of gray level intensity (focused from 140 to 165), consistent with Figure 7a,b. As the wind speed increases, the facet size decreases, and their number increases. For instance, in Figure 8c where U = 3.3 m/s, the facet size decreases rapidly while the facet quantity increases, which results in a large variant in the number of facets (with different slopes). To sum up, the shape of the histogram depends on of each facet and . It should be noted that most of the bright spots are due to the reflection of the diffused light source at random slopes in Figure 7f.

Laser Propagation
During the experiments, a high-speed camera is located on the top of the wavetank just behind the translucent receiver plate (whose height is 27 cm) to capture instantaneous images of the laser beam footprint after passing the air-water interface. The parameters considered in this study are incident angle , wind speed U, and laser beam diameter d. For each configuration, 8700 images were captured and processed to calculate the beam centroid and compare it with a free turbulent interface. The centroid coordinates were calculated based on an image moment algorithm [10,49,50]. The displacement changes in the x-direction (along-wind) or y-direction (cross-wind) from the centroid coordinate of the flat surface show how the laser beam refracted at the surface. To quantify the refraction, we converted each coordinate to a deviation angle from the flat surface's refraction angle, as shown in Figure 6. To better investigate the influence of the ocean facets on the laser propagation at different wind speeds, a collimated 0.3 cm beam was used as the input beam into the beam expander system to obtain an expanded beam of 0.8 cm. Samples of the original beam and the expanded beam are displayed in Figure 9.
The shape of the laser spots in Figure 9c,d for which U = 0.9 m/s are comparable to Figure 9a,b for the quiescent interface. The drift of the laser spot can be observed when the wind speed increases, as shown in Figure 9e,f, for which U = 1.9 m/s, and Figure 9g,h, for which U = 3.3 m/s. Based on the laser beam footprints and their deviation, the variation in the laser beam centroid location was calculated. Figure 10 shows the centroid distribution of the laser beam at θ i = 0 o at different wind speeds for two different values of d.
The centroid location is slightly changed at minimum wind speed in Figure 10 a,b and medium wind speed in Figure 10c,d due to the capillary waves on the interface. However, when the wind speed increases, the wavelength increases, and gravity waves are generated. Therefore, the centroid locations drift rapidly in Figure 10e,f compared to those resulting from the minimum and medium wind speeds. Moreover, the deviation in the along-wind direction is higher than in the cross-wind. Figure 11 presents the standard deviation of the centroid drift at different wind speeds at =

Laser Propagation
During the experiments, a high-speed camera is located on the top of the wave-tank just behind the translucent receiver plate (whose height is 27 cm) to capture instantaneous images of the laser beam footprint after passing the air-water interface. The parameters considered in this study are incident angle θ i , wind speed U, and laser beam diameter d. For each configuration, 8700 images were captured and processed to calculate the beam centroid and compare it with a free turbulent interface. The centroid coordinates were calculated based on an image moment algorithm [10,49,50]. The displacement changes in the x-direction (along-wind) or y-direction (cross-wind) from the centroid coordinate of the flat surface show how the laser beam refracted at the surface. To quantify the refraction, we converted each coordinate to a deviation angle from the flat surface's refraction angle, as shown in Figure 6. To better investigate the influence of the ocean facets on the laser propagation at different wind speeds, a collimated 0.3 cm beam was used as the input beam into the beam expander system to obtain an expanded beam of 0.8 cm. Samples of the original beam and the expanded beam are displayed in Figure 9.
The shape of the laser spots in Figure 9c,d for which U = 0.9 m/s are comparable to Figure 9a,b for the quiescent interface. The drift of the laser spot can be observed when the wind speed increases, as shown in Figure 9e,f, for which U = 1.9 m/s, and Figure 9g,h, for which U = 3.3 m/s. Based on the laser beam footprints and their deviation, the variation in the laser beam centroid location was calculated. Figure 10 shows the centroid distribution of the laser beam at θ i = 0 • at different wind speeds for two different values of d.
The centroid location is slightly changed at minimum wind speed in Figure 10a,b and medium wind speed in Figure 10c,d due to the capillary waves on the interface. However, when the wind speed increases, the wavelength increases, and gravity waves are generated. Therefore, the centroid locations drift rapidly in Figure 10e,f compared to those resulting from the minimum and medium wind speeds. Moreover, the deviation in the along-wind direction is higher than in the cross-wind. Figure 11 presents the standard deviation of the centroid drift at different wind speeds at θ i = 0 • . The standard deviation is almost equal for a beam diameter of 0.3 cm in Figure 11a and a beam diameter of 0.8 cm in Figure 11b at the minimum and medium wind speeds. On the contrary, since the beam size in Figure 9a is smaller than the facet size at the higher winds, as pointed out in Section 2.1 in Figure 7f, the standard deviation for the along-wind and the cross-wind for a beam diameter of 0.3 cm in Figure 11a is double that exhibited by the 0.8 cm beam diameter in Figure 11b. 0 . The standard deviation is almost equal for a beam diameter of 0.3 cm in Figure  11a and a beam diameter of 0.8 cm in Figure 11b at the minimum and medium wind speeds. On the contrary, since the beam size in Figure 9a is smaller than the facet size at the higher winds, as pointed out in Section 2.1 in Figure 7f, the standard deviation for the along-wind and the cross-wind for a beam diameter of 0.3 cm in Figure 11a is double that exhibited by the 0.8 cm beam diameter in Figure 11b.  0 . The standard deviation is almost equal for a beam diameter of 0.3 cm in Figure  11a and a beam diameter of 0.8 cm in Figure 11b at the minimum and medium wind speeds. On the contrary, since the beam size in Figure 9a is smaller than the facet size at the higher winds, as pointed out in Section 2.1 in Figure 7f, the standard deviation for the along-wind and the cross-wind for a beam diameter of 0.3 cm in Figure 11a is double that exhibited by the 0.8 cm beam diameter in Figure 11b.   Another way of looking at the along-wind and cross-wind displacements of the laser beam centroid due to waves is shown in Figure 12. The histograms there show, on the ordinate, the number of images (that is, the number of centroids) in which the laser beam centroid had a deviation of a value given on the abscissa. Due to the direction of the water flow, which moves parallel to the wind direction as described in Figure 5, the deviation of the displacement centroid from the quiescent interface along the wind direction outpaces the deviation in the cross-wind axis. Figure 12c,d show that the fluctuation increases as the wind speed increases. For U < 1.9 m/s, the histograms in the along-wind and the cross-wind directions nearly follow a Gaussian distribution as shown in Figure 12a Another way of looking at the along-wind and cross-wind displacements of the laser beam centroid due to waves is shown in Figure 12. The histograms there show, on the ordinate, the number of images (that is, the number of centroids) in which the laser beam centroid had a deviation of a value given on the abscissa. Due to the direction of the water flow, which moves parallel to the wind direction as described in Figure 5, the deviation of the displacement centroid from the quiescent interface along the wind direction outpaces the deviation in the cross-wind axis. Figure 12c,d show that the fluctuation increases as the wind speed increases. For U < 1.9 m/s, the histograms in the along-wind and the cross-wind directions nearly follow a Gaussian distribution as shown in Figure 12a,b. Figure 11. The standard deviation in the along-wind direction (σ x ) and cross-wind direction (σ y ) direction with 0.3 cm beam diameter (a) and 0.8 cm beam diameter (b) at different wind speeds at Figure 11. The standard deviation in the along-wind direction (σx) and cross-wind direction (σy) direction with 0.3 cm beam diameter (a) and 0.8 cm beam diameter (b) at different wind speeds at = 0 .
Another way of looking at the along-wind and cross-wind displacements of the laser beam centroid due to waves is shown in Figure 12. The histograms there show, on the ordinate, the number of images (that is, the number of centroids) in which the laser beam centroid had a deviation of a value given on the abscissa. Due to the direction of the water flow, which moves parallel to the wind direction as described in Figure 5, the deviation of the displacement centroid from the quiescent interface along the wind direction outpaces the deviation in the cross-wind axis. Figure 12c,d show that the fluctuation increases as the wind speed increases. For U < 1.9 m/s, the histograms in the along-wind and the cross-wind directions nearly follow a Gaussian distribution as shown in Figure 12a,b. Contrastingly, the deviation of the displacement centroid for maximum wind speed (U < 1.9 m/s) does not fit a Gaussian distribution. Instead, a Gaussian Kernel-smoothing distribution [51] provides the best fit, as shown in Figure 12c,d. Furthermore, the laser beam diameter impacts the deviation of the displacement centroid according to the facet size, as discussed earlier. For instance, the beam centroid with a diameter of 0.8 cm deviates from the still surface case by ±0.02 cm for U = 0.9 m/s in both the along-wind and the cross-wind directions. The deviation for maximum speed is measured as ±1 cm and ±0.5 cm in the along-wind and the cross-wind directions, respectively. For the laser beam diameter of 0.3 cm, the deviation increases. The fluctuation of the centroid measured in the along-wind direction is ±3 cm.
When the laser beam passes the air-sea interface through a single facet, its new direction and path are subject to the along-wind and cross-wind slopes belonging to that facet. The deviation of the refraction angles θ r and slope angles θ s in the along-wind axis and the cross-wind axis are calculated and obtained by a centroid shift between each sample and the original quiescent interface, as discussed in Section 2. The standard deviation of the slope angles and the refraction angles for two different beam diameters (0.3 cm and 0.8 cm) and two different angles of incidence (0 • and 32 • ) at different wind speeds are shown in Tables 2 and 3. These statistical results are in good agreement with those of a previous report [9,10] in which the deviation of both the refraction angles θ r and slope angles θ s in the along-wind direction was found to be higher than in the cross-wind direction. Based on the centroid measurements that were discussed in the previous section, the temporal power spectral density (PSD) was obtained in order to further analyze he motion of the air-sea interface. Each realization consists of 8734 time samples (frames), and these samples were divided into eight subsets, each consisting of 1024 samples. By using the fft function in MATLAB, the average of these PSDs was then computed. Temporal spectra (both along-wind and cross-wind) of the interface using the lowest and highest wind speeds and the beam diameters of 0.3 cm and 0.8 cm are shown in Figure 13.  Figure 13 shows that the relative power density of the along-wind component always higher than that of the cross-wind components, especially at lower frequenc Further, the laser beam diameter does not change the PSDs calculations qualitative For example, Figure 13a,b represent the PSD at minimum wind speed for beam dia eter d = 0.3 cm and d = 0.8 cm, respectively. Their peak energy is located approximat at the same frequency, which is f~4 Hz. Similarly, Figure 13c,d represent the PSD at m imum wind speed (3.3 m/s) for beam diameter d = 0.3 cm and d = 0.8 cm, respectively. Th peak energy is located at approximately the same frequency, which is around f~10 Hz.
The spectra shown in Figure 14 are for the interfacial slopes and approximated in facial elevation. The equations in Figure 15 (see [9,44,[52][53][54][55]) incorporate a linear appr imation to obtain the approximated 1-D elevation wave spectrum ( ).  Figure 13 shows that the relative power density of the along-wind components is always higher than that of the cross-wind components, especially at lower frequencies.
Further, the laser beam diameter does not change the PSDs calculations qualitatively. For example, Figure 13a,b represent the PSD at minimum wind speed for beam diameter d = 0.3 cm and d = 0.8 cm, respectively. Their peak energy is located approximately at the same frequency, which is f~4 Hz. Similarly, Figure 13c,d represent the PSD at maximum wind speed (3.3 m/s) for beam diameter d = 0.3 cm and d = 0.8 cm, respectively. Their peak energy is located at approximately the same frequency, which is around f~10 Hz.
The measurement duration for each experiment presented in Figure 5 was about 8 s (8700 frames of data). The interfacial slope spectrum and the corresponding interfacial elevation spectrum for our data were compared to the PM and JONSWAP spectra presented in Table 1 by applying similar parameters in Figure 5. At the lower wind speed, the peak frequency for the slope and the elevation spectra for our data are close to both models, as in Figure 14a,b, respectively, owing to the spectral peak enhancement factor for JONSWAP to improve the spectrum accuracy at the lower frequency [24] The peak frequency for our data tends to increase as the wind speed increases. Contrarily, the peak frequency of PM's model and JONSWAP's model decrease as the wind speed increases. The temporal slope spectrum can be converted approximately from the temporal elevation spectrum by (2πf) 4 g 2 . Hence, the slope spectrum falls off slower than the elevation spectrum at the higher frequencies.
The temporal and spatial frequency conversions for both the slope and elevation spectra are shown in Figure 15.
at the same frequency, which is f~4 Hz. Similarly, Figure 13c,d represent the PSD at maximum wind speed (3.3 m/s) for beam diameter d = 0.3 cm and d = 0.8 cm, respectively. Their peak energy is located at approximately the same frequency, which is around f~10 Hz.
The spectra shown in Figure 14 are for the interfacial slopes and approximated interfacial elevation. The equations in Figure 15 (see [9,44,[52][53][54][55]) incorporate a linear approximation to obtain the approximated 1-D elevation wave spectrum ( ).  The measurement duration for each experiment presented in Figure 5 was about 8 s (8700 frames of data). The interfacial slope spectrum and the corresponding interfacial elevation spectrum for our data were compared to the PM and JONSWAP spec-  The measurement duration for each experiment presented in Figure 5 was about 8 s (8700 frames of data). The interfacial slope spectrum and the corresponding interfacial elevation spectrum for our data were compared to the PM and JONSWAP spectra presented in Table 1 by applying similar parameters in Figure 5. At the lower  Figure 15. Temporal and spatial frequency conversion for both the slope and elevation spectra.

Conclusions
In this work, we conducted a set of experiments to study dynamic oceanic waves and their influence on laser propagation. We investigated the statistical behavior of irregular small-scale, air-sea interface facets under different configurations. The results show that the facet area decreases when the wind speed increases, which results in raising the slope angle of each facet. In addition, a water-air optical system was built to investigate the propagation of laser beams impacted by the turbulent interface. The results indicated that the offset deviation in the along-wind direction escalates more than in the cross-wind direction with the wind speed. As a result, the standard deviations of the refraction angle, as well as the slope angle, increase.
One limitation of this study that will be addressed in future work was the maximum wind speed of only 3.3 m/s. Improving this model by investigating the laser propagation with bigger beam diameters (larger than the facet size) at the interface with higher wind speeds will enable reliable models which apply to the corresponding larger dynamic range, especially at very lower frequency regions. Moreover, we will design and implement an optical communication link to evaluate the BER performance under those conditions in future work.
In general, the wavy air-sea interface causes deterioration of the laser communication channel. The present study emphasizes the need to investigate techniques to mitigate that effect. Active adaptive optics should be considered to reduce the impact of wavefront distortion [56]. This approach can be extended, for example, by using micro-electromechanical system (MEMS) [57] deformable mirrors.
These characteristics of the oceanic facet's behavior combined with laser propagation would facilitate future research to improve the performance of free space optical communication through the air-sea interface by tracking the optimal facet with the desired path to the receiver.