Multi-Sensor Approach to Improve Bathymetric Lidar Mapping of Semi-Arid Groundwater-Dependent Streams: Devils River, Texas

: Remote sensing technology enables detecting, acquiring, and recording certain information about objects and locations from distances relative to their geographic locations. Airborne Lidar bathymetry (ALB) is an active, non-imaging, remote sensing technology for measuring the depths of shallow and relatively transparent water bodies using light beams from an airborne platform. In this study, we acquired Lidar datasets using near-infrared and visible (green) wavelength with the Leica Airborne Hydrography AB Chiroptera-I system over the Devils River basin of southwestern Texas. Devils River is a highly groundwater-dependent stream that ﬂows 150 km from source springs to Lake Amistad on the lower Rio Grande. To improve spatially distributed stream bathymetry in aquatic habitats of species of state and federal conservation interest, we conducted supplementary water-depth observations using other remote sensing technologies integrated with the airborne Lidar datasets. Ground penetrating radar (GPR) mapped the river bottom where vegetation impeded other active sensors in attaining depth measurements. We conﬁrmed the accuracy of bathymetric Lidar datasets with a di ﬀ erential global positioning system (GPS) and compared the ﬁndings to sonar and GPR measurements. The study revealed that seamless bathymetric and geomorphic mapping of karst environments in complex settings (e.g., aquatic vegetation, entrained air bubbles, riparian zone obstructions) require the integration of a variety of terrestrial and remotely operated survey methods. We apply this approach to Devils River of Texas. However, the methods are applicable to similar streams globally.


Introduction
Remote sensing technology uses active or passive propagated signals (e.g., light pulses, electromagnetic radiation) to enable remote acquisition of information about objects.Remotely sensed data are geospatial in nature because they are collected by global positioning system (GPS) bearing instruments that allow rapid data capture over large areas [1,2].Since the late 1960s, traditional aerial photographic systems have been replaced with airborne and spaceborne sensors capable of acquiring data at different spatial, spectral, temporal, and radiometric resolutions, providing the means for a variety of analytical opportunities.
Airborne and satellite-based remote sensing approaches show great promise for characterizing aquatic ecosystems [3,4].One such approach is light detection and ranging (Lidar), an active, non-imaging, remote sensing technology that uses an artificial light source commonly falling in the optical spectrum between 0.5 and 1.6 µm.The applications and usage of Lidar technology have increased since the 1960s and cover a wide range of military, civilian, and scientific activities [5][6][7].A typical non-stationary Lidar system is composed of a laser scanner, transmitting and receiving units, a GPS, an inertial navigation system (INS), a digital storage unit, a digitizer, and an operator interface.These system components emit pulsed or continuous laser energy and with each pulse measure, the distance travelled between the transmitter, the target surface, and back to the receiver.Onboard electronics digitize the received electromagnetic (EM) energy and record the amplitude, the geographic location, and various other information.
Airborne Lidar bathymetry (ALB) is a derivative of Lidar technology, and popularly utilized for mapping relatively shallow and transparent bodies of water and shorelines from an airborne platform [8,9].A typical ALB system uses two laser transmitters: a near-infrared (NIR, 1064 nm) wavelength to detect water surfaces, and a visible (green, 515-532 nm) wavelength to measure water depths.In theory, green wavelength beams emitted from an airborne laser transmitter traverse the air-water interface and propagate in the water until reaching and illuminating the bottom.The light beams slow while traveling in the water column, and their energy (amplitude) attenuate due to refraction and scattering.If a beam is not attenuated entirely, it reflects from the bottom, and the backscatter reaches the receiver.A digitizer records the reflected amplitude in a waveform fashion, indicating the surface, bottom, and all other major interactions, including noise.Contrary to green wavelength beams, depending on the turbidity levels and surface characteristics, NIR beams are either absorbed immediately beneath the water surface or reflect off, indicating the surface location (Figure 1).Especially in transparent waters, advanced data-processing algorithms can be combined to extract elevations from both NIR and green-wavelength returns, producing a more accurate surface representation.
In recent years, researchers have completed a number of inland bathymetric surveys using ALB technology.In doing so, they have demonstrated this technology's versatility and accuracy [10][11][12][13].The term "inland" refers to specific ecosystem types, such as lakes, reservoirs, rivers, ponds, swamps, and wetlands.The term also represents ecological and hydrologic environments of diverse shapes and sizes and various physical, chemical, and optical properties [14].Especially in fluvial karst environments, the environmental conditions and aquatic properties vary greatly, presenting an accuracy challenge for airborne Lidar bathymetry [15].Therefore, it is essential to conduct in-situ measurements to quantify bathymetric expectations [16,17].With the advent of recent survey technologies, researchers now have access to autonomous unmanned surface (USV) and aerial vehicles (UAV), sophisticated GPS/INS positioning, dual and multibeam sonars, submerged spectrometers, and Lidar fluorosensors to measure depth, water transparency, bottom properties, and environmental conditions that impact inherent optical properties (IOP).
The present study addresses the data gap of stream bathymetry, which is required to parameterize physical stream-flow models to assess aquatic habitat conditions under a range of stream-flow conditions.In addition to thermal and hydrologic data collection, accurate stream bathymetry is critical to hydraulic model development.To this end, we applied a suite of complementary remote sensing approaches to characterize distributed stream depth.To map the stream bathymetry, we conducted airborne Lidar surveys and performed supplementary measurements to build a seamless map of the river bottom and immediate surroundings.
The study proposed a detailed evaluation of the pools by Dolan Falls because of their significant role in habitat conservation of the Devils River minnow.We defined the pool mesohabitat areas upstream of Dolan Falls as "upper basin" and the area downstream of the falls as "lower basin."The pool of the upper basin includes areas with dense submerged variable-leaf watermilfoil (Myriophyllum heterophyllum) and emergent American water willow (Justicia americana), which pose challenges for Lidar and sonar.Because of access issues to the watershed (98% is privately owned), in situ measurements were limited and clustered by The Nature Conservancy's Dolan Falls Preserve Area (29 • 53 ′ N, 100 • 59 ′ W). Figure 2 illustrates the entire study area and exemplifies upper and lower basins, the location of in-situ measurements (L1-11), GPS base stations, and Dolan Falls.

Materials and Methods
In this project, we generated bathymetric calculations of Devils River by integrating several survey methods.We acquired airborne Lidar datasets using a Leica Airborne Hydrography AB (AHAB) Chiroptera-I system with a maximum effective range of 500 m above the water surface.The system acquired and recorded green-wavelength data with a pulse rate of 35 kHz in a waveform format.A near-infrared scanner produced Lidar pulses at 260 kHz-optimum speed at an average survey altitude of 450 m.The Chiroptera-I has two integrated cameras: a medium-format 50-megapixel (MP) Hasselblad camera (red-green-blue, RGB) for orthophoto production, and a lower-resolution camera (video graphics array; VGA) to assist airborne missions and data processing.The Chiroptera-I is rated as a shallow (i.e., depth penetration to 15 m in relatively transparent waters) mapping system.The pulse power and energy are lower (0.062 MW and 0.12 mJ) than a typical deep-channel (i.e., depth penetration to 30 m) ALB system, such as the Leica HawkEye III (0.72 MW and 2.3 mJ).The beam divergence of each pulse is between 0.5 and 3.0 mrad (NIR, green).The Chiroptera-I uses two Palmer scanners to distribute Lidar pulses in an elliptic pattern; the maximum incident angle is ±14 • fore and aft and ±20 • toward the edges (Figure 3).Palmer scanning mirrors revolve inclined in a cycle, reflecting laser pulses to ground and projecting footprints overcast ground along spiral line [18].The advantages of a Palmer scanner include the capability to map sloped surfaces and measure the terrain beneath moderate tree canopies such as those usually found in riverine survey locations, particularly at riparian zones.As an unconventional hydrographic survey instrument, we used a ground penetrating radar (GPR) to characterize the depth and river bedrock where surface or submerged vegetation prevented other active sensors such as Lidar and sonar from accurately measuring its depth.In order to have higher chances of depth measurement and easier navigation with a boat, we completed GPR surveys in October with lesser vegetation presence due to seasonal changes.GPR surveys are common in seismic and near-surface profiling of bedrock [19,20], emitting electromagnetic pulses in the microwave band of the radio spectrum (UHF/VHF) and detecting layer boundaries or objects in the subsurface with different permittivity.EM waves are discharged and received by an antenna mounted immediately above the water surface and are recorded when the instrument is in the rover (moving) mode.Because EM waves do not travel at a constant speed through different types of material or through the water column, the accurate propagation velocity needs to be calculated based on the electrical and magnetic properties of the survey environment.The relative permittivity of the soil, rock, water, and air vary greatly.Typically, EM waves travel 30 cm/ns in a vacuum and 3.3 cm/ns in freshwater [21].

Study Site
Devils River is a semi-arid karst stream in southwestern Texas in which most discharge originates as spring flows and distributed groundwater inflows [22].The river is pristine due to its remote and rugged location, bordering primarily by large intact ranchlands and protected natural areas, which minimizes pollution from human and animal waste [23].Riparian vegetation includes stands of Mexican white oak (Quercus polymorpha) and American sycamore trees (Platanus occidentalis); upland areas are dotted with scrub juniper (Juniperus) and mesquite (Prosopis velutina).The study area is comprised of several important pool mesohabitats with a diversity of hydrologic habitats, supporting 28 species of fish [24].The protected Devils River minnow (Dionda diaboli) is classified as a federally threatened species.D. diaboli is a minnow species (family Cyprinidae) with a distribution limited to tributaries of the middle Rio Grande in Texas and Coahuila, Mexico, with adults reaching approximately 5-6 cm in length.Texas Parks and Wildlife Department (TPWD) biologists are working with conservation and research partners to develop and execute a multidisciplinary research program to collect the scientific evidence needed to inform management of D. diaboli and several other species of interest to state and federal conservation efforts.This includes data collection to support development of a robust two-dimensional hydraulic fish habitat model, enabling resource managers understand the declines in spring flow and surface water discharges that affect habitats for imperiled fish and mussel species.Findings would allow for data-driven revisions to current instream flow recommendations for Devils River.As part of this research program, researchers used an infrared camera mounted to an UAV to acquire stream surface temperature, identify spring discharge, and characterize stream and spring water mixing [25].Additionally, in-stream stage and temperature measurements were acquired to quantify spring inputs and identify seasonal trends in thermal regime affecting target species [26].

Data Acquisition
A Cessna 206 survey aircraft, owned and operated by the Texas Department of Transportation (TxDOT), spent approximately 30 h over the survey site (Figure 4A,B).The airborne survey consisted of 189 flight lines and covered an approximate area of 220 km 2 .The NIR sensor produced an average point density of 10-12 per m 2 with an approximate flight speed of 120 knots.However, we computed an average ground-sampling density of 50 points per square meter, which is highly inflated due to a number of flight lines overlapping each other, especially at river bend locations (Figure 2).To comply with quality control procedures for ALB surveys [17] and to build a continuous map of the river bottom, we acquired GPS, sonar, and GPR data at pools upstream and downstream of Dolan Falls.Ground-truthing of river water surface and bottom elevations were completed using a Trimble R8 global navigation satellite system (GNSS) receiver mounted to a variable height pole.To augment spatial coverage of the bathymetric survey, we installed a consumer-grade Garmin echoMap 74dv to a kayak and coupled it with a submerged GoPro HERO5 camera (Figure 5).The sonar encompasses a dual-band transducer with an effective range of 36 cm to 700 m in freshwater bodies, and the unit samples positional information at a frequency of 5 Hz (Figure 6).Prior to fieldwork, we tested the unit in a pool at Applied Research Laboratories (ARL), The University of Texas at Austin (30 • 23 ′ 05 ′′ N, 97 • 43 ′ 34 ′′ W).The equipment measured depths with a bias of 2.7 cm (standard deviation to the mean) at a depth of 11.8 m.Because the sonar beam diameter is much narrower than the Lidar pulse, we expected deeper measurement capability in areas with aquatic vegetation.For example, a green-wavelength Lidar pulse emitted at a range of 450 m will cover an area of 135 cm in diameter on the water surface, expanding to an area of 250 cm at 5 m below the surface.However, the sonar transducer will generate a much smaller footprint of 86 cm in diameter at the same depth, increasing the likelihood of penetration through submerged vegetation [27].For GPR surveys, we used a GSSI SIR-3000, coupled with a GSSI 200 MHz shielded monostatic antenna.The main purpose of using the GPR was to investigate the feasibility of EM waves passing through submerged aquatic vegetation and detecting the river bottom (Figure 7).The propagation of the GPR waves through different layers of materials or surfaces creates varying and distinctive resistivity in the frequency.The phase velocity (v) and attenuation constant (α) are the main parameters describing the EM wave propagation, and they both depend on the dielectric and conductivity (σ) properties of the medium.In theory, the lower the propagation (the higher the permittivity), the less the lateral dimension of the medium is resolved [28,29].Because water is highly polarizable and has a higher degree of permittivity (ε), it absorbs EM waves quicker as the antenna frequency increases.Freshwater permittivity is commonly referenced as 80, attenuation constant (α) as 0.1 dB/m and conductivity (σ) as 0.5 mS/cm [21,30].EM velocities are usually referenced in cm/ns and frequencies in MHz; therefore, assuming a 200-MHz antenna and a water temperature of 20 • Celsius, wavelength (λ) is 16.77 m and propagates (ν) at 3.3 cm/ns.For optimum results, the GPR antenna should have no interference and should be clear of signal-blocking obstructions.Therefore, we mounted the instrument to a boat with a deflated rubber bottom, and the operator did not carry any electronics onboard.We connected an external Trimble R8 GNSS receiver to the GPR and synchronized the time between the instruments for all surveys.(Figure 8).Additional in-situ measurements, such as moisture and salinity content, are suggested in the presence of mixed materials [23], but they are not considered critical in freshwater reservoirs such as the Devils Rivers watershed.

Optical Properties of Water
Light beams interact with a number of constituents while traveling through the water column.Foremost, water in its purest form exhibits a complex absorption spectrum, and a significant amount of scattering occurs due to refractive index fluctuations [31].The two fundamental IOPs of water are the absorption coefficient and the volume scattering function.Apparent optical properties (AOPs) depend on the medium and the geometric (directional) structure of the ambient light field.According to Smith and Baker [32], AOPs are particularly important when considering the penetration of radiant energy to varying depths.ALB uses green wavelengths (515-532 nm) because of the least water absorption (attenuation) (Figure 9).
In shallow waters and localized surveys, limnologists typically use conventional survey methods and refer to Secchi disk depth (Z sd ) to quantify water-column transparency.The diffuse-attenuation coefficient K d (z, λ) is an important optical property.It defines the logarithmic depth derivative of ambient irradiance E d (z, λ), which comprises photons heading in downward and upward directions [33].K d (or K PAR ) depends on the directional structure of the ambient light field and is classified as an AOP.It is determined in the photosynthetically active range of 400-700 nm.
Freshwater contains certain quantities of dissolved organic matter (DOM) and decayed vegetation that are responsible for temporal and spatial variability of optical properties.In the water column, these materials are the major cause of light and radio-wave absorption and scattering.Although there have been 21st-century advances in optical electronic instruments using tungsten light sources to measure water transparency and calculate K d [34], a popular, simple method is to use a 30-cm-diameter Secchi disk.The Secchi-disk depth (Z sd , m) represents the maximum depth at which the disk is no longer observed as it is progressively lowered into a water column.Using a Secchi disk to determine K d is simple: K d = α/Z sd .However, this measurement is generalized and may not be accurate for large and multiple pools.Lee et al. [35] defined the relationship between Z sd and K d and reported that α ranges between ~1.3 and 2.0, as observed in various studies for over 90 years [36,37].In this study, we assumed α = 1.6 and K d < 1.5 to calculate the theoretical depth measurability (D max ) of the Chiroptera-I as specified by the manufacturer [38].
A turbidimeter measures the loss of intensity in the transmitted light resulting from scattering by DOM, which vary in shape, color, and reflectance and range from <1.0 to several 10 s of micrometers in size.Researchers observe and report turbidity either in nephelometric turbidity units (NTU) or in formazin nephelometric units (FNU).In this study, we used a HACH 2100Q portable turbidimeter for observations in NTU, and the instrument calibration indicated compliance with the manufacturer recommendation (Supplementary Table S1).We sampled the river at an average depth of 1 m at 11 locations and replicated the practice three times in each location to comply with U.S. EPA 180.1 [39].We observed an average turbidity of 2.76 NTU, and the standard deviation of all readings was 0.42 NTU, which indicates suitability for Lidar bathymetry [17].

Data Analysis
Lidar bathymetry is observed by calculating the total discrete time signals (∆t) occurring between the water surface (t 1 ) and the bottom peaks (t 2 ) indicated in the backscattering signal.Variations in water transparency, depth, and other environmental conditions may limit or prevent measurement, affecting the modeling of the backscatter amplitude and the shape of the recorded waveform.Cottin et al. [41] and Eren et al. [42] investigated waveform shapes, backscatter amplitudes, and attempted to model the bottom properties.For river surveys, Allouis et al. [43] and Pan et al. [44] studied the waveform shape, modifying it via decomposition and deconvolution methods and comparing their findings to in-situ depth measurements.
Because most NIR beams get absorbed at the surface or in the immediate water column, only a few backscattering signals return to the receiver.Using a receiver with a high signal-to-noise error ratio, it is possible to detect and sort the returns and generate a representation of the water surface.The Chiroptera-I's manufacturer-supplied software, Lidar Survey Suite (LSS), outputs returns to separate classes with respective wavelength type, shape, and amplitude.Class 0 includes returns from both wavelengths, and Class 5 stores solely green wavelengths.A previous study reported a robust goodness-of-fit between these classes (R 2 = 0.99) in the transparent and deep waters of the lower Colorado River [45].Conversely, in Laguna Madre [17], which is a large, shallow lagoon in southern Texas, the goodness-of-fit between these classes declined (R 2 = 0.32) due to increased amounts of turbidity (2-10 NTU) and poor transparency (K d > 2).Therefore, variance in environmental conditions represents one of the great challenges of ALB surveys, requiring frequent and localized in-situ measurements to predict measuring-depth performance.
As for water-bottom detection, LSS outputs two separate classes of green-wavelength returns: regular (Class 7) and enhanced (Class 10).The latter output algorithm filters out weaker signals and selects the returns omitted by the Class 7 algorithm that are recorded below an optimum threshold.
In the Colorado River study, results revealed depth detection improved by an average of 0.81 m (9.2%) at less turbid conditions (0.4 NTU) and locations where Class 7 returns were extinguished.
Lidar point-cloud data were output in LAS v1.2 format, and we converted water-surface and water-bottom classes into text format (xyz) for analysis in Matlab v17.LSS requires amplitude and backscatter thresholds to be optimized for accurate bathymetric datasets, especially at survey locations where environmental conditions are dynamic.To comply with quality assurance procedures, the field staff investigated the local conditions carefully and prepared supplemental information (e.g., pictures, field notes) to assist data analysts.
Point-to-point correspondence or comparison with Lidar point-cloud datasets [46] or any other geospatial vector datasets is not assumed; therefore, we used the Delaunay triangulation algorithm [47] to create surface patches of triangulated irregular networks (TIN) (Surface 1, Figure 10).In theory, the algorithm maximizes the minimum angles of triangles in the triangulation process by defining an empty circumcircle [48].Every TIN patch represented an averaged elevation computed from the returns that formed a triangle (Surface 2, Figure 10).We compared the TIN patches to sonar, GPS, and GPR-derived elevations that registered proximity (dS1 i ) of 1.0 m and excluded the returns at slopes (α) steeper than 45 • .The vertical threshold (h) was set at 0.5 m to prevent the algorithm from picking up vegetation or other highly sloped surface returns.We performed the data analysis in ellipsoidal height elevations, which is native to GPS measurements.This is a recommended quality control procedure for a reason: the geographic positioning information acquired by most survey equipment (e.g., airborne Lidar, sonar, GPR) requires conversions of ellipsoidal heights to real-world elevations (orthometric height), and the interpolation process can carry errors in extended baselines-particularly in river and shoreline surveys.It is considered good practice to transform elevations following the completion of vector data alterations (e.g., data cleaning, reclassifications) because most Lidar and remote sensing projects require rasterized products in orthometric elevations.An average geoid height of −22.87 m was calculated using the GEOID12B model near Dolan Falls.
For GPR data analysis, it is possible to calculate the velocity (ν) of an EM wave (C em ) in an insulator by factoring the permittivity (ε w ) into the speed of light (Equation ( 1)).The wavelength (λ) in any medium equals the velocity divided by the EM frequency (Equation ( 2)).We computed the attenuation coefficient (α) using the observed water-column conductivity (σ) and the calculated permittivity (Equation ( 3)).We adjusted the EM propagation in the water column by refining the permittivity value with observed in situ water temperature (T) using Equation ( 4) as defined by Klein and Swift [49].Furthermore, we calibrated the GPR using point-source reflectors that produce "hyperbolas" in the data, which are caused by the cone shape of the antenna's radiation pattern.We fitted the hyperbolas with shapes of known velocity to the hyperbolas recorded at Devils River.Therefore, we were able to confirm wave velocities and compute the correct propagation through the water column [50].
Throughout Devils River, water temperature was recorded using the sonar transducer's imbedded thermometer, and we computed the median temperature to be 24.55 • Celsius.Using an Aqua-TROLL 200 Data Logger, we measured the median conductivity at σ = 0.44 mS/cm and computed the attenuation coefficient at α = 0.08 dB/m (Equation ( 3)).Because the river depth, discharge, and temperature vary, we expected to observe fluctuations in permittivity.

Results
Ultimately, the Devils River study attempts to (1) confirm the feasibility of using a geophysical instrument to measure the river bedrock when active remote sensing technology is not able and (2) build a seamless stream map by confirming the positional bias of a number of geospatial and geophysical datasets.Furthermore, the findings of this study are expected to contribute to the ongoing investigation of airborne Lidar bathymetry quality-control aspects applicable to karstic riverine environments with varying environmental and morphological conditions.

Water Transparency and Turbidity
The river bottom was visible at five in situ locations.Hence, we expected the Lidar beams to measure the water bottom without complications.A quadratic relationship (Figure 11) was observed, and a high degree of coefficient of determination was computed between K d and D max (R 2 = 0.99, RMSE = 0.02 m).At other locations, Z sd varied from 1.4 m to 2.9 m (Figure 12).We expected the Lidar beams to attenuate rapidly in the less transparent waters.Table 1 illustrates the Z sd and turbidity observations at all locations.AHAB claims that the Chiroptera-I can measure depths to 1.5 times the K d , indicating the D max capability [38].To confirm this specification, we computed the theoretical D max = 2.72 m by Dolan Falls (L6), where sonar readings indicated a depth of 5.77 m and Lidar beams were attenuated at a depth of 2.65 m.

Water-Surface Detection Analysis
It is challenging to detect water surfaces accurately with airborne Lidar, especially in riverine environments.In the Devils River project, we confirmed the surface representation with point-cloud datasets (Class 0 and 5).Because the Chiroptera-I emits NIR beams faster and with a narrower pulse width, fewer Class 0 returns are registered.It is expected that most NIR beams are absorbed by the immediate water surface and only a few reach back to the receiver, especially in choppier and turbid waters.
In this analysis, we used Class 5 (green wavelength only) returns as the reference points to compute the point-to-patch relationship.TIN patches were formed by Class 0 returns (NIR and green wavelength), and each patch was compared to the nearest Class 5 return by using the least-squares method.Computations revealed that a higher degree of correspondence was evident in the lower basin due to less aquatic vegetation (Table 2).Furthermore, findings indicated that Class 0 returns characterized the water-surface elevations slightly higher than Class 5 returns (Supplementary Figure S1). Figure 13
We conducted eight missions to complete airborne data acquisition.In four of these, we merged TXDR and CampTPWD base-station solutions using GrafNav application to achieve a closer baseline (mean = 40.2km) and a more precise differential solution.The application can support up to 32 base station.An average of 99.34% quality-one (Q1) resolution was achieved for all missions (Supplementary Table S2).Positional dilution of precision (PDOP), which specifies the error propagation and measurement precision of satellite geometry [51], averaged at 1.61 out of 10 scale (1-best, 10-worst).Further, we conducted a forward/reverse separation analysis to estimate the robustness of solutions and achieved RMSE of 1 cm for horizontal separation (easting and northing) and 5 cm for vertical separation.
Primarily, the CampTNC base station was used to differentially correct the static in-situ measurements.We analyzed the expected bias applicable to horizontal and vertical accuracies by computing the median standard bias and outliers in the differential solutions for 26-27 June 2019 (Supplementary Materials, Figure S2).Findings revealed a comparable horizontal and vertical bias for both days (2.5 cm and 4.1 cm).The probability plot indicates normal distribution.However, vertical solution has a higher probability to include outliers that are greater than 5 cm for both days (Supplementary Figure S3).

Lidar Bathymetry versus GPS Measurements
We measured 29 water surface elevations and 487 bottom elevations and refined these static GPS measurements differentially for higher accuracy.It was essential to measure the bottom elevations, particularly at locations with dense vegetation, to assist in building a seamless rasterized map.TIN patches were created (h L ) using all available water-bottom Lidar returns (Classes 7 and 10) and were compared to GPS-derived (h g ) surface (W S ) and bottom (W B ) elevations.
Results revealed a robust coefficient of determination (R 2 = 0.86) to GPS-derived elevations; however, only a fraction of GPS measurements (21%) were compared due to aquatic vegetation blocking the Lidar beams.Overall, GPS measured 0.03 m deeper than Lidar bathymetry (Table 3), possibly due to circumference discrepancies between the measuring rod and the Lidar beam at the water bottom.We measured the vertical bias at 0.12 m (RMSE).The regression plot illustrates an improved coefficient of determination in the higher elevations of the basin (Figure 14).In our computations, GPS-derived elevations indicated a higher probability of deeper measurement in the shallower areas of the river compared to Lidar TIN patch elevations (Supplementary Figure S4).

Lidar Bathymetry versus Sonar
We processed a total of 13,720 sonar measurements in the upper and lower basins, which came out to a mean depth of 1.26 m and 2.0 m, respectively.Sonar measured the shallowest depth as 36 cm and the deepest depth as 7.87 m at the plunge pool by Dolan Falls.The density plot shows the distribution of all sonar measurements, revealing the distinctive aquatic and geomorphologic characteristics of each basin and the significant effects caused by Dolan Falls and other environmental factors (Supplementary Figure S5).
We compared sonar measurements to Lidar-derived (W B ) TIN patch elevations using least-squares algorithm.The computations matched 7240 sonar measurements (52.77%) to TIN patches that were within the defined thresholds (refer to Section 2.4).Results indicated a robust coefficient of determination (R 2 = 0.78 and σ 2 = 0.08) due to its more uniform pool structure.Mean elevation difference computed between the datasets was slightly higher in the upper basin due to the presence of aquatic vegetation (11 cm versus 9 cm).In the lower basin, a lesser determination of coefficient was observed (R 2 = 0.72) and variance was higher (σ 2 = 0.17).Nevertheless, mean elevation difference was lower due to fewer amounts of aquatic vegetation existence (Table 4, Figures 15 and 16).

Ground Penetrating Radar
Using an inflatable boat and a GSSI SIR-3000 GPR with a 200-MHz antenna, we conducted geophysical surveys in the upper basin.We created and recorded 148 survey lines.Because GPR radio waves are continuous, we were able to extract returns at every 10 cm of forward motion and compared each measurement to the nearest Lidar-derived TIN patch.Due to dense submerged vegetation, the least-squares algorithm did not register 57% of these patches (634 samples out of 1113) with a GPR measurement (refer to Section 2.4).Preliminary regression results revealed an average depth bias of 5 cm and RMSE of 16 cm using the literature suggested fresh-water EM propagation rate of ν = 3.3 cm/ns.Signal discontinuity caused by temperature and permittivity variability was evident, particularly in the deeper areas.Of all GPR measurements taken, the deepest was 3.11 m, which confirms the depth limitations of using a 200-MHz antenna [28].To improve the findings, we adjusted the water-column permittivity (ε) with the observed mean water temperature, 24.55 • Celsius.As a result, EM waves propagated quicker in the water column (ν = 3.38-3.7 cm/ns), especially in the deeper parts of the river, and produced a better determination of coefficient (Table 5, Figures 17 and 18).
As an example, and for visual confirmation, GPR measurements at survey line #117 were output and compared to Lidar bathymetry.It was evident that Lidar returns and GPR measurements matched at locations with no vegetation or less vegetation.On a normal "B-scan" (Figure 19), the horizontal axis represents distance and the vertical scale is in two-way reflection time, which we converted to depth.As expected, submerged vegetation blocked Lidar beams in the water column, and EM waves propagated through the medium, indicating the bedrock with distinctive amplitude in the waveform.

Discussion
In this study, we attempted to confirm the positional bias of remotely sensed measurements to one another before generating a complete stream bathymetric map of Devils River by the Dolan Falls Preserve area.In inland pools, especially in fluvial karst environments, conditions can be challenging and typically have a dynamic nature.The results depicted in this study reveal that mapping with airborne Lidar technology is accurate and wholesome (98%).However, complementary surveys and in-situ measurements are required to validate bias and to complete mapping efforts of the remaining survey location (2%) in building bathymetric products.
In the upper basin by Dolan Falls, aquatic vegetation prevented airborne Lidar and sonar beams from measuring depth; therefore, we applied a novel approach and conducted GPR surveys to increase the likelihood of measuring the riverbed.GPR data acquisition and analysis was the final step in this study due to uncertainties surrounding the bias errors and completeness of data measured in such survey locations.Therefore, initial efforts (Figure 20A-C) in this study confirmed the mean vertical bias of Lidar bathymetry against post-processed GPS measurements (11 cm RMSE) and dual-beam sonar measurements (27 cm RMSE).Our final efforts validated the accuracy of GPR measurements to the TIN patches generated from Lidar returns (~17 cm RMSE).Overall, the mean vertical difference was computed less than 10 cm among comparison of all datasets.Using GPR, our research efforts primarily focused on measuring the bedrock at locations with dense aquatic vegetation.Accuracy and its quantification are essential.However, the foremost concern remains the measurement reliability caused by varying depth and permittivity.Prior to the Devils River survey, University of Texas at Austin researchers conducted a demonstration survey in Austin, Texas.Barton Springs pool has similar depths and bedrock characteristics to Devils River, allowing the researchers to investigate the suitability of using GPR throughout Devils River.Preliminary results revealed that a 200-MHz antenna is suitable for measuring depths up to 3 m, confirming GPR would be suitable for the Devils River study.A slower frequency antenna may generate deeper measurements but produce fewer details.
Other influences of EM attenuation, such as varying levels of turbulence, salinity, and bottom reflectance, may have an impact on the quality of the remotely sensed data.In the Devils River survey, the water was calm and shallow with minimal turbulence.Therefore, we did not require inertial measurements to correct positional information acquired from sonar and GPR.However, for choppier and more turbulent conditions, using an inertial measurement unit (IMU) coupled with a GPS is advisable to prevent random angle bias in depth measurements from a boat.For instance, a 3 • antenna angle, compared to the nadir angle, equals an additional 16 cm at a depth of 3 m.
Still, there are uncertainties in using active remote sensing technologies for bathymetric surveys.Because of differences in platform heights (airborne and on a boat), the sensitivity and capacity of electronics, and varying pulse widths and angles, minor discrepancies are expected.Overall, the Devils River study findings revealed that these technologies could complement each other well and the observed discrepancies are insignificant given the benefits gained from building a seamless bathymetric map. Figure 20D illustrates the combined depth map of the upper basin by Dolan Falls where aquatic vegetation interfered with Lidar and sonar measurements.

Conclusions
In Devils River, the foremost goal was to measure the river's depth at all locations with airborne Lidar, and supplement Lidar datasets using other hydrographic survey technologies as necessary.Therefore, we constructed our hypothesis based on the following assumptions: • the probability of measuring depth through aquatic vegetation with all available methods, and • the confirmation that all integrated remote sensing datasets are accurate.
We demonstrated in this article that building seamless stream maps using only Lidar datasets may not be possible, and quantifying the possible measurement bias between different survey technologies is essential.Therefore, we conducted several in-situ measurements and systematically compared remote sensing datasets to each other.We confirmed the possibility of combining datasets sourced from different remote sensors, and the feasibility of using GPR for measuring depths while aquatic vegetation is present.The quantification methods demonstrated in this study can be exhaustive and complicated.However, they are required to comply with internal quality assurance and quality control procedures.
Remote sensing hardware manufacturers provide specific training and support, and we recommend collaboration before and during mapping efforts.To avoid costly re-measurements, researchers may follow previously confirmed and established in-house quality control and quality assurance procedures.For this purpose, researchers may consider implementing the U.S. Geological Survey's recently released Lidar base specifications (version 2.1) for all products used in conjunction with their 3D Elevation Program (3DEP), which is a nationwide effort [52].Various federal and local organizations have also published their procurement standards and specified quality control procedures to be followed prior to data and product acceptance (e.g., Texas Strategic Mapping Program, StratMap).
For inland ALB surveys: passive remote sensing technologies (e.g., satellite or airborne spectral imaging), can be highly beneficial for predicting environmental conditions in which in situ measurements may be inadequate or not feasible.There is a variety of tools available.Researchers should investigate the spatial and spectral resolutions of their project requirements and available budget prior to purchasing products.A new and an exciting update to the repertoire of tools for civilian use is the no-cost data availability of NASA's Ice, Cloud, and Land Elevation Satellite 2 (ICESAT-2), which became operational in September 2018.The satellite carries a photon-counting Lidar unit with a 10-KHz scanner operating with green-wavelength energy.Recent articles have demonstrated its versatility with nearshore [53,54] and reservoir studies [55].Therefore, we believe the hydrographic remote sensing community may benefit from a study that compares the survey technologies demonstrated in this article to depths measured by ICESAT-2 for quality control purposes, especially at karst riverine environments.Theoretically, positive comparison results would reduce some of the labor-intensive in situ efforts, especially at remote locations.Additionally, due to the routine orbiting of the satellite, researchers can investigate temporal changes in project locations against the impacts of varying environmental conditions.

Figure 1 .
Figure 1.Scattering and refraction occur when light beams traverse through the air-water interface.Near infrared (NIR) beams are absorbed either at the surface or in the immediate water column, whereas some reflect back to the receiver and indicate the water surface.Distances (range) from the surface and bottom can be computed by measuring the flight duration of the pulses and adjusting with the varying speed of light in air and water.

Figure 2 .
Figure 2. Survey area and the Devils River watershed by Dolan Falls.In situ locations are marked through L1-11, and purple lines indicate the airborne survey paths.

Figure 3 .
Figure 3.A representation of the Palmer scanner and its elliptical scanning pattern of Leica AHAB Chiroptera-I sensors.

Figure 6 .
Figure 6.Surveying stream depth using a dual-beam Garmin echosounder.

Figure 9 .
Figure 9. Graph indicating the minimal absorption coefficient in the 200-800 nm spectral region, attained in pure freshwater (modified after Beć and Huck [40]).

Figure 10 .
Figure 10.Conceptual representation of Lidar, sonar, GPS, and GPR-derived elevation comparisons using Delaunay triangulation and a point-to-patch relationship.

Figure 11 .
Figure 11.Quadratic relationship between K d (m −1 ) and D max (m) in Devils River, indicating the lessening D max in parallel to the escalation in K d .
illustrates the "noisy" nature of Class 5 returns compared to Class 0 output (Standard deviation = 10 cm versus 3 cm) in the lower basin, where NIR beams reflected off and penetrated the water column up to 0.19 m of median surface representation.

Table 2 .Figure 13 .
Figure 13.Representation of Class 0 (NIR + green) versus Class 5 (green wavelength only) returns, indicating the water-surface detection challenges in inland pools.Note the "noisy" characteristic of Class 5 returns; indicating the median surface is 9 cm lower than Class 0 returns.

Figure 15 .
Figure 15.Regression of sonar-derived elevations versus Lidar-derived TIN patch (W B ) elevations in the upper basin.Aquatic vegetation is evident at the surface and in deeper parts of the water column.

Figure 16 .
Figure 16.Regression of sonar-derived elevations to Lidar-derived TIN patch (W B ) elevations in the lower basin.Sonar measured deeper than Lidar due to depth by the Dolan Falls plunge pool area.

Figure 17 .
Figure 17.GPR measurements compared to Lidar TIN patches using the least-squares method.

Figure 18 .
Figure 18.Solid linear-fit line revealing an improved regression between Lidar and GPR-derived depths (3.3 cm/ns versus 3.7 cm/ns) using the adjusted EM propagation.

Figure 20 .
Figure 20.Maps showing upper basin depths measured using (A) GPS, (B) sonar, and (C) GPR, and (D), a final depth map with all measurements combined.
Figure S1: Distribution of Class 0 versus Class 5 for surface representation.The upper basin registered a greater sample range and worse deviation, Figure S2: Median horizontal and vertical accuracy bias present in differential correction, CampTNC base station (26-27 June 2019), Figure S3: Probability of positional bias in differential solution applicable to localized measurements, Figure S4: GPS-derived elevations indicate higher possibility of deeper measurement in the shallower areas of the river compared to Lidar-derived TIN elevations, Figure S5: Density plot indicating the distribution of sonar-derived depths in the upper and lower basins surrounding Dolan Falls.Author Contributions: K.S. supervised the project as co-principal investigator, performed fieldwork, analyzed datasets, and drafted the article.A.R.A. operated the Chiroptera-I ALB system, assisted with fieldwork, and contributed to sections of the article.L.C. performed fieldwork, operated the GPR, and processed GPR datasets.B.D.W. conceived the project as principal investigator, received research funding, contributed to sections of the article, and provided insight.S.R. initiated the project at TPWD and provided research funding.All authors have read and agreed to the published version of the manuscript.Funding: This study was funded by the Texas Parks and Wildlife Department (TPWD) via the U.S. Fish and Wildlife Service and the Bureau of Economic Geology, John A. and Katherine G. Jackson School of Geosciences, The University of Texas at Austin (Contract No: 507663).

Table 1 .
Secchi-disk and turbidity measurements at Devils River.Areas with dense aquatic vegetation were excluded.VB indicates a visible river bottom.

Table 3 .
Lidar-derived (h L ) water-bottom (W B ) patch elevations compared to GPS measurements (h g ).
σ σ Figure 14.Goodness-of-fit between Lidar-derived patch elevations to GPS measurements in the upper basin.

Table 4 .
Lidar-derived TIN patch elevations versus sonar measurements for the upper and lower basins by Dolan Falls.

Table 5 .
GPR measurements compared to Lidar TIN patches, indicating the mean and maximum depth discrepancies with initial and final EM velocity coefficients.