The Effect of Surface Waves on Airborne Lidar Bathymetry (ALB) Measurement Uncertainties

Airborne Lidar Bathymetry (ALB) provides a rapid means of data collection that provides seamless digital elevation maps across land and water. However, environmental factors such as water surface induce significant uncertainty in the ALB measurements. In this study, the effect of water surface on the ALB measurements is characterized both theoretically and empirically. Theoretical analysis includes Monte Carlo ray-tracing simulations that evaluate different environmental and hardware conditions such as wind speed, laser beam footprint diameter and off-nadir angle that are typically observed in ALB survey conditions. The empirical study includes development of an optical detector array to measure and analyze the refraction angle of the laser beam under a variety of environmental and hardware conditions. The results suggest that the refraction angle deviations (2σ) in the along-wind direction vary between 3–5◦ when variations in wind speed, laser beam footprint size and the laser beam incidence angle are taken into account.


Introduction
Airborne Lidar Bathymetry (ALB) is a remote sensing technology that uses pulsed green lasers to measure coastal water depths [1]. The ability to collect rapid data to generate digital elevation maps makes the ALB an important tool in mapping, surveying and hydrographic applications [2]. The laser beams that are emitted from an aircraft interact with the water surface, water column and seafloor and return to the receiver in the aircraft. These interactions are recorded as time-series signals, that is, waveforms. Bathymetric measurements are calculated by taking into account the time-of-flight difference between the water surface (surface return) and the seafloor (bottom return). However, environmental factors such as water surface, water column and seafloor induce uncertainty in the ALB measurements [3,4]. Among these factors, water surface plays a primary role on the measurement uncertainty as it significantly alters the laser's ray-path direction and geometry as the laser pulse refracts into the water column [1].
In typical ALB surveys, the water surface (air-water interface) is treated as a horizontal plane where the laser beam refracts into the water column at an angle depending on the off-nadir beam angle. In nature, however, wind-driven surface waves modify the incidence angle of the laser beam, that is, the angle between the laser beam optical axis and the water surface normal. As a result, the assumption of a horizontal plane is no longer valid and the refracted ray-path geometry of the laser beam underwater will depend on the water surface conditions. The surface waves vary in wavelength and amplitude from the millimeter scale (capillary waves-dispersed by surface tension) up to tens of meters (gravity waves-dispersed by the force of gravity). Depending on the characteristics of One of the main uses of ALB surveys is to update charts for marine navigation. The quantification of measurement uncertainty in ALB surveys is critical for charting as the surveyor is required to report the quality of the survey through uncertainty measurement. According to the International Hydrographic Office (IHO) S-44 standards publication [10], horizontal and vertical total propagated uncertainty (TPU) components of the ALB survey should be provided in addition to the bathymetry. Based on a NOAA forensic analysis between ALB and multi-beam echosounder surveys, it was considered that ALB surveys conducted by hydrographers are typically considered to be order 1b [11].
Previous studies investigated the contribution of the air-water interface on the laser pulse propagation in various aspects. These include the contribution of the air-water interface on the laser beam energy, target detection capability and pulse geometry during its travel in the environmental medium. It is shown that the air-water interface results in variations in downwelling and upwelling directions and laser beam energy losses [2,12,13]. Increasing wind speeds also impedes target detection capability that is important for navigation application and bottom analysis. It has been shown in [14] that higher wind speeds decrease the probability of detecting a 1 m 3 cube. The surface waves have also been demonstrated to increase the variance in the ray-path length for laser beams with smaller beam footprints [15]. The effect of sea state on the accuracy of ALB measurements has been discussed in various studies. It has been suggested that at strong wind conditions for ALB surveys with wind velocities of 10-12 m/s, the maximum horizontal and vertical position errors correspond to 5-10% and 1-2% of the depth, respectively [16]. The effect of sea state and the beam divergence angle on the laser beam coordinate accuracy has been demonstrated in previous studies. A water surface consisting of small wavelets (capillary and capillary-gravity waves) can cause considerable displacements to an ALB's rays in the horizontal and vertical axes [17]. A ray-tracing model was used to estimate the effects of capillary and gravity waves on the ray-path geometry [18]. This was followed by an empirical analysis based on a single sea state condition. Although previous studies discussed the effect of sea-state on the ALB measurements, we have identified two gaps that we attempt to address in this study: (1) The water surface has not been decoupled from other environmental contributions (such as from water column and seafloor) in order to calculate a TPU error budget for the ALB; and (2) extended empirical validations, needed to validate the model results, are missing.
In this study, the effect of water surface on the ALB measurements is investigated by decoupling it from water column and seafloor. Specifically, the effect of capillary-gravity waves on the laser beam path is evaluated using theoretical and empirical approaches. In the theoretical approach, Monte Carlo ray-tracing model is developed to simulate the laser beam geometry changes under varying water-surface conditions. The variation in the laser beam path is assessed through modeling water surface with wavelengths that are smaller than the transmitted laser beam footprint diameter. The empirical approach includes the use of an optical detector array that is built to intersect the laser beam's light field and image the spatial distribution of the laser beam beneath the water surface. Refraction angle variations of the laser beam are assessed from the optical detector array output using image processing algorithms. The empirical measurements are conducted in well-controlled environment at the University of New Hampshire's (UNH) Ocean Engineering facilities. Based on the empirical measurements, the total horizontal uncertainty (THU) and total vertical uncertainty values (TVU) are calculated. The results of the study provide baseline for water surface contribution to the overall ALB TPU budget.

Water Surface Model
The geometry of the air-water interface is mainly determined by the wind and surface tension conditions and is classified into two main wave types: capillary waves (i.e., waves with wavelengths smaller than 0.02 m and restored by surface tension) and gravity waves (i.e., larger swells with wavelengths larger than 0.02 m that are restored by the gravity) [19]. In this study, capillary-gravity waves are taken into account as the surface wave condition that is most likely to be observed during an ALB survey. To understand the contribution of the capillary-gravity water surface waves on the ALB measurement uncertainty, the air-water interface is modeled based on Fourier transform methods [20]. The water surface generation procedure begins with choosing a 1D wave spectrum that describes the wave regime, that is, whether the wave is capillary wave, capillary-gravity wave or gravity wave based on wavelength. The wave spectrum chosen for this study is the Apel spectrum because it describes the capillary-gravity wave regime [21]. Apel wave spectrum is a modified version of the Joint North Sea Wave Project (JONSWAP) spectrum [22] that includes an improved prediction of capillary-gravity waves. Another motivation for using the Apel wave spectrum in the resulting simulated water surface models is that sea state conditions generated from short fetch winds (less than 1 km) in waters shallower than 100 m provide a close approximation to the ALB survey conditions. The 1D Apel wave spectrum is given as [23]: where A is the spectral constant adjustment term to fit significant wave height data chosen such that A = 0.00195, k is the spatial frequency, L PM is the Pierson-Moskowitz shape spectrum defined as: k p is the wavenumber of the spectral peak that depends on the inverse wave age and the wind speed as: where k 0 is the wavenumber at nominal windspeed U 10 measured at a height 10 m above the water surface and g is the gravitational acceleration. The inverse wave age, Ω, is related to the dimensional fetch, x, measured in meters and defined based on field measurements as [22]: J P is the JONSWAP peak enhancement term defined as [22]: with γ depending on Ω as: and Γ is defined as where σ is the adjustment term for the wave age: Experimental data obtained from the wind-wave tanks were used to tune the capillary-gravity wave regime [24]. Accordingly, the high-pass and resonance behavior of the spectrum were modeled. R ro is the high-pass quadratic filter with roll-off wavenumber of 100 rad/m: and R res is the resonance gravity peak defined as: where constants a = 0.8, k res = 400, k w = 450 and sech is the hyperbolic secant. s is the curvature level of an assumed secondary gravity-peak at 750 rad/m: V dis is the dissipation mechanism observed in the experimental spectrum and modeled as: with the dissipation frequency k dis = 6283 rad/m. I D defines the integration of the vector spectrum over the azimuth direction: where φ s = 0.28 + 10 k p k 13 −0.5 Through the modeling steps, the 1D Apel wave spectrum in the spatial frequency domain is converted into 3D water surface elevation [20]. The environmental information, such as wind speed and fetch, are used as inputs to the Apel spectrum model to provide 3D water surface elevation models that approximate the water surface conditions during an ALB survey [25]. Because it is inherently complex to measure a 2D wave spectrum, the approach to determining the wave direction is by using directional spreading functions [26]. The spreading function estimates the spread of possible wave propagation directions relative to direction of the wind rather than assuming that all the waves propagate along the same direction. A cosine-2S wave spreading function is used to describe the directional component of the wave spectrum as: where f is the frequency, ψ is the angle away from the mean wind direction angle, ψ o , γ is the Gamma function and s p is the dimensionless spreading parameter that varies from 1-15. The s p value was chosen as 2, as high frequency waves related to capillary-gravity waves have spreading parameter of 2 [27]. The full 2D directional wave spectrum, is calculated by taking the product of 1D Apel wave spectrum with the cosine-2S spreading function as follows: The Apel wave spectrum, cosine-2S spreading function, and the resulting directional spectrum created with wind speed of 5 m/s at a fetch length of 5.5 m are provided in Figure 2. The final step for simulating the water surface is the creation of normally distributed random Hermitian amplitudes. The resulting Hermitian amplitudes are multiplied by the 2D directional wave spectrum to ensure that each entity in the 2D directional wave spectrum has an inverse Fourier transform. Taking the inverse Fourier transform of the spectrum results in 3D water surface elevations. The generated water surface is triangulated using Delaunay triangulation because a triangular surface provides facets that are used to calculate the refraction angle of the laser ray ( Figure 3). Then, the water surface facets along with the water surface normal are calculated and are used in the ray-tracing model. It should also be noted that the described model does not take into account hydrodynamic forces such as currents and waves formed by water displacement. The reader is referred to Reference [20] for more details on the water surface modeling procedure.

Ray-Tracing Model
Numerical modeling is used to identify key parameters for empirical analysis as a limited amount of ALB survey configurations can be conducted in laboratory conditions due to the logistical constraints. After a realistic water surface elevation model is constructed over the spatial domain, a total of 10,000 simulated light rays (composing a laser beam) are traced through the water surface to calculate the ray distribution at a depth of 0.25 m below the water surface. The depth of the detector array is selected to minimize the effects of water clarity, that is, scattering and absorption, on the ray-path geometry. The simulation parameters are selected to match the empirical measurements conducted at UNH Ocean Engineering facilities. In the ray-tracing simulation model, the laser source is assumed to be a point source and is located at 4 m height with mean water surface elevation assumed to be at 0 m. The laser beam divergence angle is kept at 60 mrad to obtain the empirical laser beam diameter footprint size of 0.25 m on the water surface. The laser rays emitted from the source are then intersected with the water surface facets and the direction vector for the refracted vector, s 2 , is determined by using the vector form of the Snell's law as follows [28]: where s 1 is the direction vector of the incident ray, N surf is the normal vector of the water facet, η a is the refractive index of air (η 1 = 1.00), η w is the refractive index of water (η w = 1.33) ( Figure 4). The environmental input parameters used in the simulations are wind speed, U 10 , and fetch, x, whereas the hardware related input parameters are laser beam incidence angle and laser beam footprint diameter on the water surface ( Table 1). The wind speed used in the simulations range from 2 m/s to 5.25 m/s. Fetch is varied from 3.5 m to 8.5 m at 0.5 m increments. The laser beam incidence angle is varied from 0 • to 20 • and the laser beam footprint diameter is varied from 0.25 m to 4 m which are typically used in both topo-bathymetric lidar and bathymetric lidar systems. A total of 500 water surface models is generated in the Monte Carlo ray tracing simulations for each case. For each water surface realization, a sub-surface laser beam centroid location is calculated by tracing the laser rays that intersect at z = 0.25 m. The reference laser beam centroid is calculated for a laser beam that passes through a smooth, flat surface (still water surface). The centroid shift between the disturbed water surface and a reference still water surface conditions indicate refraction angle deviation caused by the capillary-gravity surface waves Figure S1. The centroid shift obtained in distance measures is converted to a refraction angle deviation value-∆θ-from the still water assumption in the along-wind and cross-wind axes. It is possible to calculate the angular offset in the laser beam refraction angle in the along-wind and cross-wind directions for each image using the following equation and the depth of the optical detector array (i.e., 0.25 m): where ∆θ aw is the laser beam refraction angle in the along-wind direction, Y i is the x-axis centroid location for the given image, and Y still is the centroid for the still water reference case. Similarly, for cross-wind direction, the laser beam refraction angle is calculated as follows: where ∆θ cw is the laser beam refraction angle in cross-wind direction, X i is the y-axis centroid location for the given image, and X still is the centroid for the still water reference case.

Data Acquisition and Image Processing
Two types of experiments are conducted to empirically measure the effect of water surface on the laser beam's ray-path geometry. The goal of the first type of experiments is to provide wind speed and wave spectrum conditions similar to the water surface models described in Section 2.1. Wind over the water tank surface is generated using an industrial fan and wind speed is measured with an anemometer (Digital Tools Mini anemometer). To characterize the spectrum of the fan-generated capillary-gravity waves, water surface elevation is measured with a capacitive wave staff (Ocean Sensor Systems OSSI-010-002) at a sampling rate of 30 Hz.
The goal of the second type of experiments is to characterize the interaction of the laser beam with different water surface conditions using the optical detector array designed at the University of New Hampshire Center for Coastal and Ocean Mapping [29,30]. The optical detector array consists of a 6 × 6 photodiode array (Thorlabs PD1A) arranged in a planar square grid with dimensions of 0.25 m × 0.25 m. The spacing between the photodiodes is kept at 0.05 m to avoid optical cross talk between the elements [30]. Each photodiode is mounted in a clear acrylic waterproof housing. To increase the dynamic range, the photodiodes are connected via SubMiniature-A (SMA) cables to a reverse bias circuitry, which supplies a reverse bias voltage of 5 Volts. The detector array intersects the refracted laser beam and samples the light field at a sampling rate of 20 Hz using 10-bit (0-1023 digital counts) analog-to-digital converters. The output of the detector array is a 6 × 6 image. The laser beam centroid is calculated by using image moment invariants algorithm with sub-pixel accuracy.
Image moment invariant algorithm is an image processing technique that reveals information regarding the symmetry of a given image [31]. The moments of images algorithm are calculations of the weighted average of the image intensities as: where M pq is the image moment invariant with p and q are index values (p = 1, 2, 3 and q = 1, 2, 3), x o and y o are the index locations of the pixel of maximum intensity, I is the pixel intensity, S is the summation of pixel intensities, that is, S = ∑ i, j I i,j , and i and j are the number of pixels along the x-and y-directions. The image centroids, x and y, which are used to compare the simulated and empirical laser beam centroids, are calculated as: In this study, the image moment invariant algorithm is used to characterize the digitized laser beam image to known geometrical conditions, such as laser beam angle of incidence and beam diameter, through a calibration procedure. The calculated image centroids provide a measurement of the laser beam centroid shift as a function of the input parameters. Similar to the Monte Carlo ray-tracing simulations, variation of the image centroid locations is compared to a reference still water surface condition in order to measure the water surface refraction deviation when the water surface is disturbed by wind.

Experimental Setup
The empirical experiments were conducted in University of New Hampshire's Ocean Engineering facilities using a wave tank (36.5 m long, 3.6 m wide and 2.4 m deep) (Figure 5a). Wind speed and water surface elevation data were collected at distances from 3.5 m to 8.5 m away from the industrial fan at 1 m increments. The positioning for data collection was provided using an actuated tow carriage with positioning accuracy of 0.01 m. A 5-min window sample was collected for statistical analysis at each configuration in order to capture a significant number of wave crests (more than 100 samples). The optical detector array was mounted on 80/20 aluminum frames to the tow carriage with an L-frame configuration and oriented horizontally to face upwards towards the water surface (Figure 5a). The modularity of the setup allowed geometrical changes in the configuration of the mount required to collect data at changing incidence angles. The detector array was submerged to a depth of 0.25 m below the water surface to minimize the scattering and absorption effects as described in Section 2.2 (Figure 5b). A 5 mW continuous wave laser pointer was mounted to the tow carriage at a height of 4 m. A plano-concave 75 mm focal length lens was used to expand the laser beam footprint diameter to 0.25 m on the water surface, which was limited mainly by the dimensions of the detector array. The laser beam incidence angle was varied from 0 • to 20 • , specifically, (0 • , 10 • , 15 • and 20 • ). The laser beam image data obtained from the optical detector array were collected for 20 s at 20 Hz, that is, 400 frames were recorded at each configuration. The output of the optical detector array at two consecutive time steps in the same experimental configuration (θ = 15 • and distance from the fan is 3.5 m) under disturbed water surface conditions is demonstrated in Figure 6a,b. The optical detector array allows the changes in the laser centroid location as well as the focusing and defocusing effect caused by the surface waves to be observed [2] Video S2.

Wave Spectrum Results
The wave-staff measurements provided time-series of wave elevation data. Using Fourier Transform techniques, wave-staff measurements were processed to obtain the experimental wave spectrum. Table 2 summarizes the key parameters obtained from the capacitive wave staff measurements that include peak frequency, wavelength and significant wave height, calculated as average of the highest one-third of the waves, for each dataset.
The results show that the peak frequency decreases linearly (R 2 = 0.93) as the distance from the fan increases within the tested range. This decrease in the peak frequency indicates that the wavelengths increase linearly with increasing fetch. Also, the frequency content of the waves remains relatively constant for all distances within 6-10 Hz range. This constant frequency content shows that the capillary waves with small amplitude and relatively high frequency content are present in all cases. The shoulder effect observed between 10-12 Hz region, that is, slight increase in power spectrum amplitude, is attributed to the noise in the signal as it approaches the limit of the sensors measurement range of waves less than 0.002 m in height. The significant wave height also increases until 5.5 m and then decreases to 0.009 m and stays constant. In order to relate the wave conditions in laboratory setting to real-world capillary-gravity wave conditions, a fully-developed Apel spectrum model is fitted to the experimental wave spectrum results. The fully developed Apel spectrum was chosen under the assumption that in real-world survey conditions, capillary-gravity waves almost always will be fully developed, that is, the waves are not limited by fetch. As demonstrated in [32], small waves in the capillary-gravity wave regimes are fully developed by a fetch of 30 m, a distance easily achieved even in near shore applications. Accordingly, an Apel wave spectrum model with a fetch of 30 m was fitted to the experimental spectrum ( Figure 7). These models were then used in the air-water interface modeling in the Monte Carlo ray tracing simulations and the results were directly compared to the empirical measurements.

Monte Carlo Simulation Results
Using Monte Carlo simulations, it is possible to quantify the effect of wind speed, laser beam incidence angle and laser beam footprint diameter on the variation of a laser beam centroid in the along-wind (axis parallel to the direction of wind) and cross-wind (axis perpendicular to the direction of wind) directions ( Table 3). The ray tracing results conducted with a laser beam diameter of 0.25 m indicate that the beam refraction angle standard deviation shows a linear pattern with wind speed in the along-wind direction (the lowest R 2 = 0.95 at 0 • incidence angle). In addition, the along-wind direction refraction angle deviation is positively correlated with the laser beam incidence angles. The maximum average incidence angle variation is observed for the 20 • incidence angle results with 3.9 • . These results indicate that ALB systems with circular or elliptical scanning patterns and off-nadir angles and that range between 15 • to 20 • are expected to have higher refraction error uncertainty in the along-wind direction than ALB systems with rectilinear scanning patterns, the off-nadir angles of which are smaller. In the cross-wind direction, no correlation was observed between the refraction angle standard deviation and the wind speed. The laser beam refraction angle standard deviations showed a constant pattern with respect to wind speed. Additional sets of simulations were conducted in order to evaluate the effect of varying laser beam footprint diameter on the variation of the refraction angle. The simulation results conducted at a 15 • incidence angle show that the larger diameter beams result in less variation in the refraction angle in the along-wind direction (Figure 8). These results verify the findings presented in [18] which state that the wave refraction error is higher for smaller diameter beam footprints than larger diameter beam footprints. Additional simulations were conducted for a range of laser beam footprint diameters that are generated by commercial ALB systems (0.25 m, 0.5 m. 1 m, 2 m and 4 m) with wind speeds of 2 m/s, 2.5 m/s, 3 m/s, 3.5 m/s, 4 m/s and 5.25 m/s and the beam incidence angles ranging from 0 • to 20 • at 5 • intervals. The results were averaged to examine the effect of beam diameter on the along-wind and cross-wind refraction angle deviation (Figure 9). The results in Figure 9 show that the variation in refraction angle for both the along-wind and cross-wind decreases with increasing beam diameters. The variation of the refraction angle in the along-wind direction is higher than in the cross-wind direction. It is observed that for ALB systems operating with a laser beam diameter of 2 m or larger, refraction angle variations in the along-wind and cross-wind directions stay relatively constant.

Empirical Results
An optical detector array was used to collect laser beam footprint images under a variety of water surface and laser beam geometry conditions, including varying wind speed, laser beam incidence angle and distances between the laser beam footprint and the fan. The refraction angle variation was assessed by using image centroid calculations described in Section 3.2. A total of 400 image frames were collected and statistically analyzed to estimate the refraction angle. Data sets were collected for laser beam incidence angles of 0 • , 10 • , 15 • and 20 • . The variation in the laser beam centroid location was calculated as a result of interaction with the water surface. The results are compared to the reference centroid calculations conducted using a still water surface condition ( Figure 10). The results in Figure 10 indicate that the standard deviation (2σ) of the beam center for the wind-driven water surface cases is around ±0.015 m at a depth of 0.25 m, which is approximately an order of magnitude larger than the standard deviation observed in the still water case (the maximum error is~0.00086 m at 2σ for along and cross-wind directions). The results also showed that the laser beam center fluctuates slightly more in the along-wind direction than the cross-wind direction due to the nature of the wave spreading in the wave tank. The observed shifts in the laser beam center of concentration is converted to refraction angles. The standard deviation values of the refraction angles obtained from the empirical measurements are demonstrated in Table 4.  The empirical refraction angle calculations showed that the along-wind refraction angle standard deviation is linearly correlated with the wind speed (the lowest R 2 = 0.84 is observed for a 15 • incidence angle, considering the standard deviation result for the case of a 0 • incidence angle and a 2.5 m/s wind speed as an outlier). However, no correlation with wind speed was observed in the cross-wind directions. Comparison between the results obtained from the Monte Carlo simulations and the empirical measurements shows similar trends and values for the refraction angle standard deviation, especially in the along-wind direction. The refraction angle standard deviation is linearly dependent on the wind speeds in both simulations and empirical measurements. The maximum average standard deviation value of 3.9 • is observed for 20 • incidence angle, which is approximately the same with the simulation results (3.8 • in simulations). Similar to the simulation results, the refraction angle standard deviation in cross-wind direction does not show a correlation to wind speed or incidence angle.
The main difference between the simulation and empirical results were observed in the numerical values of the refraction angle standard deviation in the cross-wind direction. The cross-wind standard deviation values obtained from the empirical measurements is on average approximately 1.5 • higher than in Monte Carlo simulations.

Discussion
The results obtained in this study can be used to estimate the THU and TVU components at various depths. According to IHO, the maximum allowable THU and TVU in Order 1b standards are defined as: In order to verify the NOAA forensic analysis that ALB surveys can be considered IHO S-44 order 1b [11], the THU and TVU values are estimated for ALB survey conditions with a wind speed of 5.25 m/s at 20 • incidence angle (Table 5). A Monte Carlo simulation was conducted to calculate depth that range from 1 m to 10 m and the horizontal displacement of the laser beam using 10,000 realizations. The empirical the along-wind and cross-wind refraction angles demonstrated in Table 4 were used to calculate the TVU and THU in the Monte Carlo simulations. Accordingly, the THU values were calculated as follows: The extended ray paths due to the refraction angle standard deviation are calculated as: where r ext is the extended ray path due to the refraction angle standard deviation. Then, the TVU values were calculated by multiplying the laser ray path with cosine of the incidence angle as: The results in Table 5 show the estimated horizontal and vertical uncertainty results estimated from the Monte Carlo simulations as well as the maximum allowable THU and TVU as defined by IHO Order 1-b. The results indicate that both estimated THU and TVU values are within the IHO Order-1b standards. It is also observed that the maximum allowable TVU is more stringent than the maximum allowable THU. For example, the calculated TVU value (2σ) at 10 m depth is~0.27 m whereas the IHO Order-1a specification is~0.52 m. These findings are critical in a sense that it demonstrates the water surface contribution to the TPU budget as a function of depth. In addition, the remaining uncertainty tolerance is quantified for other uncertainty mechanisms that are not included in this paper, for example, trajectory, scanning angle, water column scattering and seafloor reflection.
The results presented in this paper can also be generalized with respect to wind speed, laser beam incidence angle and laser beam footprint size on the water surface. The refraction angle variation in the along-wind direction is linearly correlated with the wind speeds (i.e., ∆θ aw = a a U + b a ), whereas the refraction angle variation in the cross-wind direction can be approximated as a constant (∆θ cw = b c ). In addition, the refraction angle variation in both across-wind and the along-wind directions increase with larger laser beam incidence angle. Overall, the refraction angle standard deviations (2σ) in the along-wind and cross-wind directions with respect to wind speed and beam incidence angle is modeled as: ∆θ aw = (0.003θ + 0.4193)w + 0.00276θ + 1.6208 (27) ∆θ cw = (−0.0009θ + 0.0063)w − 0.00281θ + 2.3 (28) where w is the wind speed. The results obtained can also be extended to systems that use different beam footprint diameters by using a correction factor. Based on the results presented in Figure 10, the correction factors for the along-wind and cross-wind refraction angle standard deviation is given as: For laser beam foot print diameters that range between 0.25 < D ≤ 4. Although the results for the along-wind directions agreed with the simulation and empirical results, the cross-wind refraction angle standard deviation values are higher in the empirical measurements than in the simulation results. In addition, in the Monte Carlo simulations, the along-wind refraction angle standard deviation values are higher than the cross-wind refraction angle standard deviation values (Table 3) when compared to empirical results (Table 4). This was attributed to the differences in the data resulting from the simulations conducted with a set of approximations and the laboratory environment. In the empirical measurements, the waves might have reflected off the wall and changed direction until they intersected the laser beam footprint, resulting in a different wave spreading pattern. These off-directional waves might have contributed in the larger cross-wind refraction angles observed in the empirical data than in the Monte Carlo ray tracing simulations. In the Monte Carlo simulations, the surface waves follow a more uniform pattern than in the empirical measurements, which results in along-wind refraction angle standard deviation values that are higher than the cross-wind refraction angle standard deviation values. Another factor that could result in the differences in the simulation and empirical results is the method of centroid calculation in both methods. In the simulation results demonstrated in Table 3, there is a decreasing trend in the cross-wind refraction angle standard deviation values with incidence angles whereas there is no such trend in the empirical results (Table 4). Since the laser beam centroid calculation relies on the laser's ray-path geometry and energy, a potential difference in the results may occur if surface wave models in the Monte Carlo simulations do not exactly represent the empirical waves.
The sub-pixel accuracy algorithm used in the image processing methods enhanced the error resolution in the empirical measurements. Using the optical detector array, the centroid estimation error was measured to be 0.00086 m (2σ) at z = 0.25 m for the along-wind and cross-wind directions as shown in Figure 10a. This translates into refraction angle measurement error of~0.2 • (2σ) in the along-wind and cross-wind directions. Although the spacing between the photodiodes in the optical detector array is configured to be 0.05 m, sub-pixel accuracy algorithm allows to calculate the positioning errors at a higher resolution.
Although there is good agreement between the simulated and the empirical measurements, the results presented in this study should be evaluated as preliminary due to the limited number of wind speeds used in the simulation and experiments. In the empirical measurements, the main limitation was the fan power that provided the wind output to generate the capillary-gravity waves. However, in survey conditions higher wind speed values are likely to occur. At wind speeds greater than 5.25 m/s, the linear trend observed between wind speed and the along-wind refraction angle standard deviation may not be valid. Instead, a nonlinear trend may be observed. Therefore, care should be taken in extrapolating the results provided in this study to wind speeds that are larger than 5.25 m/s.
Another important contribution of the study is regarding the effect of ALB hardware on the measurement uncertainty. It is observed that larger laser beam diameters result in less refraction angle variation due to the laser beam interaction with the small wind waves. The decrease in refraction error for larger laser beam footprints is attributed to an averaging effect, that is, the total number of wave crests within the footprint is large enough that the majority of the laser rays is still refracted along the direction predicted by the still water assumption. This means that for topo-bathymetric lidar systems that use smaller footprints than conventional bathymetric lidar systems, the measurement uncertainties can be greater. Another significant observation is found for the ALB systems that use circular or elliptical scanning patterns whose off-nadir angles range between 15-20 • . The refraction angle variation at this range is found to be greater in circular or elliptical scanning patterns than the ALB systems that use rectilinear scanning patterns with smaller off-nadir angle.
Although this study provided a method to decouple the uncertainties caused by the air-water interface, further investigation is needed to produce a comprehensive ALB uncertainty model that includes the effects of uncertainties caused by both hardware and environmental uncertainties. These uncertainties include the aircraft trajectory, scanning unit, atmospheric effects, scattering in the water column and turbidity as well as seafloor contributions. The results obtained from the variation of these uncertainty parameters can be used in a look-up table to report horizontal and vertical TPU on per-pulse basis.

Conclusions
In this study, the contribution of the air-water interface on the ALB measurement uncertainty is assessed. Monte Carlo ray tracing simulations along with empirical measurements were used to measure refraction angle variation as a result of laser beam interaction with the air-water interface. The results show that the wind speed is a dominant parameter that affects the refraction angle standard deviation. In addition, the wind direction has an impact on the wave refraction angle variation. The variation in the refraction angle in the direction parallel to the wind (the along-wind) is greater than in the orthogonal direction (cross-wind). Good agreement between the Monte Carlo simulations and empirical results was found. In both approaches, the refraction angle standard deviations (2σ) in the along-wind direction vary between 3-5 • when the variations in wind speed, laser beam footprint size and the laser beam incidence angle are taken into account.
The results also suggest that different hardware configurations used by different lidar manufacturers may result in different uncertainty measurements. ALB systems with circular or elliptical scanning patterns that have off-nadir angles ranging between 15 • and 20 • are expected to have higher refraction error uncertainty in the along-wind direction than ALB systems with rectilinear scanning patterns, the off-nadir angles of which are smaller. Also, smaller beam footprint systems typically employed by the topo-bathymetric lidar systems tend to have higher refraction angle variation when the dominant wave mechanism is capillary-gravity waves. The results of this study can be used to estimate water surface contribution to the overall TPU budget and provide a baseline for look-up tables that can be used in bathymetric lidar TPU reporting.