Shallow Water Depth Estimation of Inland Wetlands Using Landsat 8 Satellite Images

Abstract

Water depth affects many aspects of wetland ecology, hydrology, and biogeochemistry. However, acquiring water depth data is often difficult due to inadequate monitoring or insufficient funds. Satellite-derived bathymetry (SBD) data provides cost-effective and rapid estimates of the water depth across large areas. However, the applicability and performance of these techniques for inland wetlands have not been thoroughly evaluated. Here, a time series of bathymetry data for inland wetlands in West Kentucky and Tennessee were derived from Landsat 8 images using two widely used empirical models, Stumpf and a modified Lyzenga model and three machine learning models, Random Forest, Support Vector regression, and k-Nearest Neighbor. We processed satellite images using Google Earth Engine and compared the performance of water depth estimation among the different models. The performance assessment at validation sites resulted in an RMSE in the range of 0.18–0.47 m and R² in the range of 0.71–0.83 across all models for depths <3.5 m, while in depths >3.5 m, an RMSE = 1.43–1.78 m and R² = 0.57–0.65 was obtained. Overall, the empirical models marginally outperformed the machine learning models, although statistical tests indicated the results from all the models were not significantly different. Testing of the models beyond the domain of the training and validation data suggested the potential for model transferability to other regions with similar hydrologic and environmental characteristics.

1. Introduction

Understanding water depth dynamics is critically important across all basic and applied wetland sciences. For example, variations in water levels affect total dissolved oxygen, nutrient, and sediment levels, vegetation composition and structure, and habitat quality for wetland-dependent wildlife, and can alter ecosystem service rates [1,2,3,4,5,6]. Water depth also affects water storage capacity and overall water quality [2]. For example, deeper wetlands can store more water which regulates flooding and retains water during drought [7,8]. Conversely, shallow wetland depths often fluctuate seasonally, promoting higher wetland vegetation productivity and enhancing carbon sequestration, effectively serving as carbon sinks [9,10]. Water depth is also critical for wetland management because it can be manipulated to achieve specific management goals. For example, increasing water depth can provide optimal foraging conditions for some waterfowl species [11] or be used to control invasive aquatic plants [12]. Likewise, decreasing water levels has been used to control mosquito populations [13] or encourage the germination of certain plant species [14]. Surface water extent and wetland depth dynamics are key factors mediating ecosystem services for wildlife and people, which generates trillions of dollars to global economies and enables systems’ resilience to climate and land use changes [15,16,17,18]. Therefore, real-time estimation of wetland extent and water depths is essential for informed decision-making; however, quantifying variation in surface water extent and water depth is challenging because monitoring systems lack needed spatiotemporal resolution and field sampling is often inefficient and cost-prohibitive.

Traditionally, surface and sub-surface wetland water levels are monitored using installed gauges, piezometers, submersible pressure transducers, and manual graduated rods [19]. These wetland water level monitors are often inadequate for specific monitoring or management goals because their spatial and temporal resolution is limited due to logistical and budgetary constraints, which is insufficient for long-term analysis [20]. Additionally, since most water level monitoring is point-based, it is difficult to extrapolate water levels for all locations within a wetland without a dense network of gauges. Spatial bathymetry can be acquired from depth sounding boats but these methods are mostly applied to deep water coastal areas, lakes, and reservoirs [1]. However, accurate boat-based measurements are limited to boat tracks leaving out very shallow waters (0–1 m) which require specific compact echo sounders [21,22]. Airborne Light Detection And Ranging (LiDAR) can achieve high spatial and vertical resolutions of ≤1 m and ≤15 cm, respectively, but LiDAR methods are costly which prohibits its use for large-scale surveys [21].

Multispectral satellite imagery has shown promise in estimating surface water extent and depth and may serve as an alternative or a complementary data stream to traditional gauge stations, echosounders, and LIDAR-derived depth estimates. In comparison, surface water and depth estimates from multispectral imagery are inexpensive, encompass vast spatial extents with repeatable coverage, and can often reveal water and depth at locales where manual measurements are impractical or inaccessible [21,22,23]. Despite the clearly identified benefits of satellite-derived bathymetry (SDB), most studies that have employed methods or algorithms that infer water depth from satellite imagery have been carried out in deep water coastal environments with relatively clear water [21,22,23,24,25,26,27]. To the best of our knowledge, no SDB method has been tested on inland shallow wetland environments that often contain turbid water mixed with vegetation and macroalgae. Furthermore, most studies have been limited to a few scenes of multispectral imagery without considering the effect of temporal variations of spectral reflectance on water depth estimation.

To overcome the aforementioned challenges, we developed and tested a method for measuring water coverage and depth in inland shallow-water wetlands through time (i.e., multiple scenes) using remotely sensed data. Specifically, we used the cloud-based Google Earth Engine (GEE) computing platform to access analysis-ready Landsat 8 satellite imagery [28] and generated bi-weekly to monthly SDB estimates using 5 empirical and machine-learning methods to study wetlands located along tributaries of the Mississippi River floodplain. We compared bathymetry estimates among empirical and machine-learning algorithms based on their accuracy, derived from water level loggers and manual water depth measurements, and interpreted measures of accuracy and precision based on biology and physical properties of satellite reflectance in shallow, turbid wetlands. The transferability of the methods we developed was tested on inland wetlands in the Florida Everglades. To the best of our knowledge, this is the first study to estimate SDB for inland shallow wetlands and provides insight into the best methods that will enable managers to use satellite imagery to estimate inland water depths across vast geographies.

2. Materials and Methods

2.1. Research Area

The study sites consist of wetlands located within the Mississippi Alluvial Valley and associated tributaries in the western part of Kentucky and Tennessee, USA. Figure 1 presents the location map and wetland sampling locations used for model training, validation, and testing. The wetlands are part of the United States Department of Agriculture’s (USDA) Wetland Reserve Program (WRP). The Natural Resources Conservation Service (NRCS) evaluates the performance of WRP restoration practices on private lands enrolled in the Program [29]. The hydrological restoration practices included the construction of shallow water levees, water control structures, levee breaks, ditch plugging, and tree planting. The wetlands are often seasonally flooded with water present for extended periods, especially in the growing season (April–October).

2.2. Water Depth Data and Image Processing

Water level data were collected from 10 wetlands in Kentucky for depth estimation model training and validation. Each of the 10 sites had an installed HOBO MX2001 water level logger that measured water levels every 15 min. Data collection occurred from August 2018 through October 2020. Seven sites were used in training the SDB algorithms and validated at three sites. The data were filtered based on the available satellite images. The average water level within 1 h before and after satellite image acquisition time was computed for each site and used to train and validate the SDB algorithms. The total number of training and validation datasets coinciding with Landsat 8 images were 266 and 126, respectively. The training datasets are used to fit and tune the parameters of the models. Two wetlands located in Tennessee were used as testing sites. Pressure transducer water level loggers have inherent errors due to measurement errors of pressure, temperature and relative humidity. Field observations also revealed errors associated with rapid drying and wetting of the pressure sensor. Depth estimation errors using in situ measurements collected by HOBO loggers were not evaluated due to resource limitations; however, these errors have been shown to be minimal [4].

Table 1. List of training, validation, and testing sites with the number of satellite images per site and geographic locations.

Site	Latitude	Longitude	No. of Images
Training site 1	36.6264	−89.1430	43
Training site 2	36.6278	−88.9615	40
Training site 3	36.6298	−88.9492	43
Training site 4	36.6151	−89.0290	43
Training site 5	36.7714	−88.9160	37
Training site 6	36.7001	−88.8016	19
Training site 7	36.9348	−88.9368	41
Validation site 1	36.6069	−89.1179	43
Validation site 2	36.6117	−89.0331	43
Validation site 3	36.7751	−88.9108	40
Testing site 1	36.2369	−89.2059	41
Testing site 2	36.1664	−89.4009	23
P36	25.5272	−80.7956	35
G-3272	25.6649	−80.5389	35
NP201	25.7166	−80.7194	35

shows the coordinates of the water level loggers.

Additional datasets were collected in Tennessee (Figure 1) by in-field manual measurement of water depths using a tape measure between December 2019 and February 2020 and repeated between November 2020 and February 2021 and November 2021 and February 2022. The manual measurements recorded in Tennessee were used as a testing dataset to ascertain the robustness of the depth models to a region outside the training and calibration data. The testing points were averaged for locations where more than one measurement was taken within a satellite image pixel. The in-field manual measurement test sites are not included in the list in Table 1.

To assess the transferability of the models in different regions outside the Mississippi alluvial plains, the models were tested at three randomly selected freshwater sites in Everglades National Park, Florida. The sites; P36, G-3272, and NP201 together with their location data can be found in Table 1. The water level data for the Everglades were downloaded from the Everglades Depth Estimation Network (EDEN) database through the Explore and View EDEN (EVE) web application (http://sofia.usgs.gov/eden). Water-level gauges in the EDEN network are operated and maintained by multiple agencies that are responsible for the accuracy of the measured data.

For the period that in situ water level data were available, cloud-free Landsat-8 images were retrieved from the Google Earth Engine (GEE) platform using PyGEE-SWToolbox [30], a Python-Google Earth Engine toolbox for surface water analysis. The toolbox automates the pre-processing and download of satellite images from the GEE platform. Pre-processing includes cloud masking, mosaicing of different satellite scenes for the same date of acquisition, and clipping of images to the study area. The pre-processing workflow was applied on atmospherically corrected Landsat 8 Level 2, Collection 1 surface reflectance data (GEE ImageCollection ID: LANDSAT/LC08/C01/T1_L2). The Landsat dataset was chosen because of the availability of data for long-term analysis compared to the recent Sentinel-2 dataset.

All the satellite images were sampled using a cloud threshold of 20%. Clouds were masked using the QA_Pixel band. At every training and testing location, we extract the reflectance values from the blue, red, green, near-infrared (NIR), shortwave infrared 1 (SWIR₁), and shortwave infrared 2 (SWIR₂) bands from all images in the image collection for the study period. The Normalized Difference Water Index (NDWI) [31] was used to extract water pixels from each image. Pixels with NDWI values greater than zero were considered water pixels. The reflectance values are log-transformed and extracted at pixels coinciding with in situ measurement locations. The satellite data acquisition and processing methodology is summarized in Figure 2. The log-transformed reflectance values become the input or predictor variables in the depth estimation models while the measured water depth is the target variable. The target variable was estimated at each pixel and evaluated at pixels coinciding with the location of the validation and testing data.

2.3. Satellite-Derived Bathymetry (SDB) Models and Modifications

All SDB algorithms can be classified into three categories: (1) empirical methods, which are the most commonly applied methods that use in situ water depth observations to correlate subsurface water reflectance to observed water depth; (2) physics-based algorithms, which do not require in situ depth measurements but use radiative transfer theories to estimate water depth; (3) machine learning methods [22]. The following sections describe the empirical and machine learning methods used in this study. The physics-based algorithms were not used in this study because their accuracy is dependent on the radiometric accuracy of input satellite images, making them highly sensitive to atmospheric correction [21].

2.3.1. Emprical Modesl and Modifications

The most commonly used empirical method for SDB estimation is the Lyzenga model Lyzenga [32] because it is efficient and simple [21]. The Lyzenga model employs multiple linear regression to determine the relationship between observed depths and log-transformed reflectance in the blue and green bands. The Lyzenga model is defined in Equation (1).

$$Z=a_0+a_i{X}_i+a_j{X}_j$$

(1)

where Z is the depth, $a_0$ is the intercept or regression constant, $a_i$, and $a_j$ are the regression coefficients determined using known water depths and corresponding reflectances through multiple linear regression, and $X_i$ and $X_j$ are the log-transformed reflectances of the blue and green bands, respectively. The independent variables are defined by Equation (2).

$$X=ln[R_w\left(\lambda \right)-R_{\infty }\left(\lambda \right)]$$

(2)

where $R_w$ is the reflectance in band $\lambda$ and $R_{\infty }$ is the average optically-deep water reflectance in the same band.

The primary limitation of the Lyzenga model is that it cannot be used to estimate depth in some conditions where X is undefined when $R_{\infty }>R_w$ is negative. However, this often can occur in inland wetlands because they are characterized by suspended solids (i.e., turbidity) near the water’s surface, resulting in greater deep-water reflectance values. Grasses and macroalgae also absorb more light which leads to lower reflectance values compared to deep-water [25]. Both ecological processes common in inland wetlands often produce negative values when subtracting deep-water reflectance from the band reflectance which makes X (Equation (2)) undefined. The model also assumes a uniform bottom albedo which prevents the model from accommodating variable substrate [21]. To address the limitations of Lyzenga’s model in areas of turbid waters, Bramante et al. [33] adapted Lyzenga’s model for the application of turbid waters by ignoring deep-water values when calculating X (Equation (3)). The Bramante version of the Lyzenga model is used in this study.

$$X=ln\left[n{R}_w\left(\lambda \right)\right]$$

(3)

where n is a fixed constant ($n=1000$) that ensures a positive log transformation and a linear response as the retrievable water depth increases.

Stumpf et al. [34] developed another empirical model for SDB estimation using the ratio of log-transformed green and blue bands to account for turbidity and improve upon the Lyzenga Lyzenga [32] model. The Stumpf model assumes a linear decrease in log-transformed ratios of high to low absorption bands and water depth [25]. The model is expressed mathematically as Equation (4).

$$Z=m_1{\displaystyle \frac{ln\left[n{R}_w\left({\lambda }_i\right)\right]}{ln\left[n{R}_w\left({\lambda }_j\right)\right]}}-m_0$$

(4)

where $m_1$ and $m_0$ are the slope and intercept determined by linear regression between the ratio and known depths, and ${\lambda }_i$ and ${\lambda }_j$ are the green and blue band spectral reflectances, respectively. The constant n has the same definition as in Equation (3).

The modified Lyzenga (Bramante et al. [33]) and the Stumpf et al. [34] models have largely been evaluated in coastal locations where turbidity levels are anticipated to be lower than inland wetland environments. Wetlands have water that is often mixed with vegetation and suspended sediments. When water contains suspended sediment, the blue band is inadequate for bathymetry estimation [35] because soil has the highest reflection (lowest absorption) at the SWIR₁ wavelengths range. With the assumption that turbidity in inland wetland waters is largely due to soil particles, we used the SWIR₁ instead of the Blue Band in the empirical models to retain the premise that the transformed ratio between the high and low absorption band and water depth should decrease linearly.

Figure 3b shows that the reflectance of SWIR₁ decreases linearly with depth with a stronger Spearman rank correlation, r = −0.59 compared to the Blue band with r = 0.15 (Figure 3a). This indicates that the SWIR₁ band is more suitable for depth estimation compared to the Blue band. The replacement of the Blue band with the SWIR₁ as used in the Stumpf band ratio model showed a stronger correlation as seen in Figure 3c,d.

To support the proposal of replacing the blue band with the SWIR₁ band, a test was conducted to select the best band combination that yielded the highest accuracy. An iterative band combination process was carried out with the six selected bands following the method by Amrari et al. [22] which yielded 15 band combinations. The band combinations were evaluated using the in situ training and validation data with mean absolute error (MAE) metric to select the optimal combination with the minimum error. The relationship between the satellite reflectances with the field-measured depths was estimated using linear regression models developed with the Scikit-learn [36] Python package to determine the parameters of the Lyzenga and Stumpf models.

2.3.2. Machine Learning Models

The Lyzenga and Stumpf models, and our modifications to them, assume a linear relationship between the water depth and the band reflectance or band ratio. In reality, the relationship may not always be linear, hence the need to explore nonlinear functions to map bathymetry. This study used machine learning models to capture these potential nonlinear relationships. The machine learning techniques considered were Random Forest (RF) regression [37], Support Vector Regression (SVR) [38], and the k-Nearest Neighbor (kNN) regression [39]. Tonion et al. [40] used both RF and SVR on Landsat, Sentinel-2, and Planetscope datasets to estimate SDB in a coastal environment and achieved reasonably accurate results.

Random Forest is a widely used supervised classification and regression machine learning algorithm that uses an ensemble of decision trees that are aggregated for optimal classification or regression predictions [41]. RF classification is for categorizing data into groups, while RF regression is for predicting continuous numerical outcomes. In the case of RF regression which is used in this study, the final prediction is the average of all the numerical values predicted by the decision tree. Each decision tree is trained on a different subset of data which helps the model to generalize better to new data. RFs are widely adopted and powerful for prediction tasks; unlike many other ML algorithms, RFs allow users to evaluate the relative importance of predictor variables, can handle many inputs and overcome multi-collinearity, and are robust against overfitting [41]. The general architecture of an RF model is depicted in Figure 4.

Support Vector Regression is based on the original Support Vector Machine (SVM) method, but instead of finding the best line/plane or a decision boundary that separates classes with the widest margin for classification, it constructs a decision boundary that lies close to as many training data points as possible that has the smallest error while allowing some flexibility to avoid overfitting [42]. The SVR algorithm maps the input space into a high-dimensional feature space using a kernel function [43] that computes the relationships in the higher-dimensional space without actually transforming the data. In the higher-dimensional space, SVR finds a linear relationship in this space that corresponds to a non-linear relationship in the original space. The radial basis kernel function is used in this study. The concept of SVR is illustrated in Figure 5 where ${\xi }_i$ and ${{\xi }_i}^{\ast }$ are slack variables that penalize training errors over error tolerance $\epsilon$. The parameter $\epsilon$ controls the width of the insensitive zone which influences the number of support vectors and thus the overall generalization ability of an SVR model [43].

The k-NN regression model is based on the k-nearest neighbors classification algorithm where an object is assigned a class based on the common class of the closest k number of training examples. The distances of nearest neighbors are often computed using the Euclidean or Manhattan distance metrics. For regression problems, the dependent variable is the mean values of its k-nearest neighbors. The average value is computed based on uniform weights where all values from all the neighbors are weighted equally or based on distance weighting where close neighbors will have a greater influence than neighbors that are further away [39].

The machine learning models were also developed using the Scikit-learn Python package. The log-transformed optimum band combination was used as predictor variables in training the machine learning models. Fine-tuning of the parameters that control the learning process of each machine learning model was conducted using a standard grid search with 10-fold cross-validation. Prior to calibrating the best hyperparameters, the optimum band combination selection for the machine learning models was carried out with the default hyperparameters for each machine learning model. Deep machine learning models were not considered in this study due to the small size of the training data.

2.4. Assessing the Effect of Turbidity on Depth Estimation

Turbid waters can affect the radiance leaving the water surface and thus, affect estimated water depth. The turbidity can increase the radiances in the visible and near-infrared spectrums, resulting in excessive depths in shallow areas and underestimating depths in deeper areas [24]. We used the Normalized Difference Turbidity Index (NDTI) [44] to quantify suspended solids and added NDTI as a predictor of water depth. The NDTI is defined by Equation (5). Positive NDTI values and low negative values closer to zero indicate turbid waters, unmasked clouds, or sun glint, while the most negative values indicate the best atmospheric and clear water conditions [45].

$$NDTI={\displaystyle \frac{{\lambda }_{Red}-{\lambda }_{Green}}{{\lambda }_{Red}+{\lambda }_{Green}}}$$

(5)

where ${\lambda }_{Red}$ and ${\lambda }_{Green}$ are the reflectance values in the Red and Green bands, respectively.

The NDTI values at each site for every satellite image were computed using Equation (5). The effect of turbidity was assessed in two ways; (1) by adding the turbidity index as a predictor variable in the depth estimation models and (2) by developing the models without the turbidity index as a first step and then combining the model results with the turbidity index in a linear regression model for final depth estimation.

2.5. Model Performance Assessment

We calculated the mean absolute error (MAE), the root mean squared error (RMSE), and the Nash–Sutcliffe efficiency (NSE) to quantify accuracy and precision for each water depth estimation method. The standard deviation (SD), Pearson’s product-moment correlation coefficient (r), and the centered pattern root mean square difference (CRMSD) are summarized in Taylor diagrams [46]. The CRMSD is the mean-removed RMSE [47]. Taylor diagrams help to visualize these statistics for all the models and the observed data in one graphical presentation. The Taylor diagrams are plotted using the SkillMetrics Python package [48]. The statistics used are expressed in Equations (6)–(10).

Differences among depth estimation methods were formally evaluated using Levene’s test with an a priori alpha = 0.05 [49]. The test was used to verify the assumption of equal variance for linear models and to determine if a significant difference exists across the models. Levene’s test was chosen because the model results were positively skewed.

$$MAE={\displaystyle \frac{1}{N}}\sum _{i=1}^N|O_i-P_i|$$

(6)

$$RMSE=\sqrt{{\displaystyle \frac{1}{N}}\sum _{i=1}^N{(O_i-P_i)}^2}$$

(7)

$$NSE=1-\left[{\displaystyle \frac{{\sum }_{i=1}^N{\left(O_i-P_i\right)}^2}{{\sum }_{i=1}^N{\left(O_i-\overline{O}\right)}^2}}\right]$$

(8)

$$r={\displaystyle \frac{{\sum }_{i=1}^n\left(O_i-\overline{O}\right)\left(P_i-\overline{P}\right)}{\sqrt{{\sum }_{i=1}^n{\left(O_i-\overline{O}\right)}^2}\sqrt{{\sum }_{i=1}^n{\left(P_i-\overline{P}\right)}^2}}}$$

(9)

$$CRMSD=\sqrt{{\displaystyle \frac{1}{N}}\sum _{i=1}^N{\left[(O_i-\overline{O})-(P_i-\overline{P})\right]}^2}$$

(10)

where $O_i$ is the observed parameter with a mean denoted by $\overline{O}$; $P_i$ is the predicted parameter with a mean denoted by $\overline{P}$, n as the number of instances.

3. Results and Discussion

3.1. Optimum Band Combination

The search for the optimum band combination for each model was carried out using all the training and validation datasets. The results of the optimum band combination assessment as shown in

Table 2. Results of band combination assessment to select optimum band combination based on MAE. Bold text indicates the best-performing model for each band combination.

Band Combination	Mean Absolute Error (MAE) (m)
	Stumpf	Bramante	RF	SVR	KNN
Blue/Red	0.948	0.934	1.042	0.893	0.949
Blue/Green	0.897	0.926	0.971	0.893	0.898
Blue/NIR	0.650	0.631	0.662	0.890	0.597
Blue/SWIR1	0.548	0.515	0.673	0.887	0.602
Blue/SWIR2	0.534	0.526	0.648	0.886	0.582
Red/Green	0.906	0.909	1.018	0.893	0.888
Red/NIR	0.679	0.655	0.662	0.889	0.624
Red/SWIR1	0.566	0.525	0.699	0.885	0.623
Red/SWIR2	0.547	0.523	0.680	0.884	0.597
Green/NIR	0.605	0.611	0.680	0.890	0.601
Green/SWIR1	0.529	0.510	0.679	0.886	0.594
Green/SWIR2	0.546	0.544	0.666	0.885	0.588
NIR/SWIR1	0.829	0.552	0.896	0.892	0.795
NIR/SWIR2	0.893	0.579	0.934	0.893	0.863
SWIR1/SWIR2	0.885	0.545	1.016	0.893	0.896

did not show a consistent optimum band combination across all models but the Green/SWIR₁ combination was dominant, particularly for empirical models that assumed linear reflectance patterns. Further analysis was carried out using the optimum band combination for each model and compared to using the dominant Green/SWIR₁ for all models to assess if the result from the optimum band combination for a particular model was significantly different from the Green/SWIR₁ at alpha = 0.05. However, no band combinations showed significant improvement in model performance so the Green/SWIR₁ band combination was adopted across all models for validation and testing.

3.2. Calibration Results

For the modified Stumpf model, an ordinary least squares regression was fitted to evaluate the relationship between the Green/SWIR₁ band ratio and the water depth training data. The overall regression was statistically significant with $R^2=0.62$ and p-value < 0.000 using the Landsat dataset. The final Stumpf model is expressed in Equation (11). A similar analysis was carried out for the Bramante (Lyzenga) model using multiple linear regression which resulted in an $R^2=0.61$ and p-value < 0.000. The final Bramante model is expressed in Equation (12). The results of the hyperparameter tuning of the RF, SVR, and KNN machine learning models are shown in

Table 3. Results of machine learning model hyperparameter tuning. These are the values of the parameters that result in the optimal performance of the models.

Model	Hyperparameter	Value
RF	No. of estimators	40
SVR	Epsilon ($\epsilon$)	0.01
	Regularization constant (C)	1
KNN	number of neighbors, k	21

.

$$Z=7.37\times {\displaystyle \frac{ln\left[n{R}_w\left({\lambda }_{green}\right)\right]}{ln\left[n{R}_w\left({\lambda }_{swi{r}_1}\right)\right]}}-6.415$$

(11)

$$Z=4.365+(0.404\times X_{green})-(0.652\times X_{swi{r}_1})$$

(12)

3.3. SDB Estimates

3.3.1. Validation Site 1

All the models performed poorly with RMSE ranging between 1.76 and 1.93 m across models. Model results were segregated into depth ranges which revealed RMSE between 0.19 and 0.29 m for depths <1 m. Depth range of 1 to 3 m also had RMSE between 0.55 and 1.30 m while depths >3 m had RMSE between 2.92 and 3.21. This indicated that the performance of the models deteriorated with increasing depth. The poor performance worsened at depths greater than 3.0 m as seen in Figure 6a. In general, empirical models outperformed machine learning algorithms based on RMSE and MAEs; the best-performing model was the Bramante algorithm (RMSE = 1.76). However, the maximum NSE across all models was 0.26, indicating poor predictive performance. In particular, we attribute the lower performance of data-hungry machine learning methods to sparse training data (

Table 4. Performance of SDB models at validation site 1 and the effect of turbidity on models. The first row for each model presents model performance without the turbidity index and the second row shows when the model result was combined with the turbidity index in a linear regression model for final depth estimation.

Model	NSE	RMSE (m)	MAE (m)	${\mathit{R}}^2$
Stumpf	0.242	1.776	1.154	0.641
Stumpf_NDTI	0.237	1.782	1.155	0.639
Bramante	0.259	1.756	1.107	0.649
Bramante_NDTI	0.237	1.782	1.155	0.639
RF	0.251	1.766	1.158	0.476
RF_NDTI	0.414	1.562	1.086	0.474
SVR	0.056	1.982	1.336	0.567
SVR_NDTI	0.491	1.456	0.969	0.565
KNN	0.105	1.930	1.255	0.574
KNN_NDTI	0.508	1.431	0.969	0.574

).

The residuals for all the models at validation site 1 are positively skewed (Figure 7a) indicating models underestimated observed depths except for the RF model. A Taylor diagram [46] (Figure 7b) is presented to visualize and compare statistical information of SDB models and how they compare with the target variable (water depth) based on their correlation, RMSE, and standard deviation. The radial distance (x and y axes) represents the standard deviation while the angles are the correlation coefficient. The distance from the reference point (observed; black dot) represented by the purple semicircles is the measure of the centered RMSE [47,50]. Good-performing SDB models will be closer to the reference point. Figure 7b shows all the models clustered away from the observed point which indicates the non-similar between the SDB estimates and observed water depths. All the model results had a standard deviation of less than 1.0 m compared to the observed data which had a standard deviation of about 2.1 m except for the RF model results which had a standard deviation of about 1.3 m. This is similar to the distribution of the residual errors seen in Figure 7a. All the performance metrics had a very small range across models (Table 4). Levene’s test indicated model results had equal variances and were not statistically different (test-statistic = 1.346 and p-value = 0.254). This means that the performances of the models are not that different despite the minor differences in the error metrics.

The turbidity index (NDTI) values obtained at validation site 1 across all sampling dates ranged between −0.44 and 0.17 with a mean of −0.07 which indicates turbid waters. The deep water samples (depth ≥ 2.0 m) had a mean NDTI of −0.03 indicating higher turbidity compared to the shallow water samples (depth < 2.0) with a mean NDTI of −0.12. The NDTI as a predictor variable did not affect any of the models’ performance. The two-step strategy of combining model results with the NDTI for final prediction demonstrated considerable improvement in prediction for machine learning models but not empirical models. The empirical models were designed to correct for water optical properties [51]; therefore, accounting for turbidity with NDTI was redundant and did not affect their performance. Conversely, at $\alpha =0.05$, the KNN showed statistically significant improvements by introducing NDTI (test-statistic = 7.716, p-value = 0.007) while SVR and RF performance did not improve significantly (test-statistic = 2.935, p-value = 0.090 and test-statistic = 0.145, p-value = 0.704, respectively). Specifically, the SVR model’s MAE improved from 1.33 to 0.97 m, whereas the MAE of the KNN and RF models improved from 1.26 to 0.97 m and 1.16 to 1.09 m, respectively, (See Table 4).

3.3.2. Validation Site 2

SDB performance at validation site 2 was much better compared to the validation site (Figure 6b). The residuals were within −0.2 m and 0.4 m with an overall pattern of underestimation (Figure 7c). The MAE obtained ranged between 0.13 and 0.19 m while the RMSE ranged between 0.18 and 0.25 m. Empirical models again outperformed the machine learning models in terms of the RMSE. Overall, the Stumpf model was the best-performing model (MAE = 0.13, RMSE = 0.18, $R^2$ = 0.80, and NSE = 0.75;

Table 5. Performance of SDB models at validation site 2 and the effect of turbidity on models. The first row for each model presents model performance without the turbidity index and the second row shows when model result was combined with the turbidity index in a linear regression model for final depth estimation.

Model	NSE	RMSE (m)	MAE (m)	${\mathit{R}}^2$
Stumpf	0.750	0.177	0.131	0.803
Stumpf_NDTI	0.747	0.178	0.133	0.803
Bramante	0.697	0.195	0.154	0.834
Bramante_NDTI	0.747	0.178	0.133	0.803
RF	0.516	0.247	0.190	0.741
RF_NDTI	0.119	0.333	0.216	0.731
SVR	0.619	0.219	0.170	0.785
SVR_NDTI	0.271	0.303	0.192	0.769
KNN	0.683	0.200	0.146	0.761
KNN_NDTI	−0.403	0.420	0.238	0.759

).

Despite different errors, we detected no statistical differences in performance among models (test-statistic = 0.669, p-value = 0.614). The superior performance of all models at validation site 2 compared to validation site 1 can be attributed to the generally shallow water depth in validation site 2. Deep water attenuates the electromagnetic signals, thereby affecting the received bottom reflectance hampering the predictive ability of water depth estimation across all models. The maximum measured water depth at validation site 1 is 6.95 m, while at validation site 2 it is 1.42 m. Validation site 1 also performed well for depths less than 2.0 m.

The assessment of turbidity at validation site 2 also indicated the waters were turbid with NDTI between −0.29 and 0.13 with a mean of −0.05. Unlike validation site 1, the inclusion of the turbidity index deteriorated all SDB model performances (Table 5). For very shallow waters, the SDB model performance is largely controlled by bottom reflectance rather than the water column properties; therefore, the turbidity index at validation site 2 was an unnecessary model parameter that worsened rather than improved performance.

3.3.3. Validation Site 3

Water depth estimates at validation site 3 had similar performance across models as in validation site 2. All the models correctly estimated the peak water depth (3.42 m) which was observed on 6 March 2019, except for the KNN model (Figure 6c). Overall, all models slightly overestimated water depths by an average of 0.27 m. Levene’s test did not show any significant differences in the results from all the models (test-statistic = 0.044, p-value = 0.996) although the KNN was the superior model based on RMSE (

Table 6. Performance of SDB models at validation site 3 and the effect of turbidity on models. The first row for each model presents model performance without the turbidity index and the second row shows when the model result was combined with the turbidity index in a linear regression model for final depth estimation.

Model	NSE	RMSE (m)	MAE (m)	${\mathit{R}}^2$
Stumpf	0.732	0.408	0.277	0.765
Stumpf_NDTI	0.728	0.412	0.281	0.763
Bramante	0.761	0.386	0.246	0.787
Bramante_NDTI	0.728	0.412	0.281	0.763
RF	0.653	0.465	0.321	0.714
RF_NDTI	0.313	0.654	0.469	0.708
SVR	0.721	0.417	0.242	0.725
SVR_NDTI	0.072	0.760	0.457	0.721
KNN	0.769	0.379	0.254	0.772
KNN_NDTI	0.237	0.690	0.427	0.769

) and balanced residual errors with few outliers (Figure 7e). Figure 7f shows that all the outputs from the models had standard deviations closer to the observed data except for the KNN which had a lower standard deviation. The waters are validation site 3 were also turbid with a mean NDTI of −0.03. Including the NDTI turbidity index deteriorated the performance of all SBD models. This result is similar to that of validation site 2 due to shallower water depths in sites 2 and 3, compared to deeper water in site 1.

3.3.4. Testing

A test was conducted to evaluate the robustness of trained SDB models in estimating depths at test locations outside the domain of the training and validation datasets (i.e., model transferability; Figure 1). The first test was carried out on the manually measured water depths in flooded corn fields where all the models performed poorly. Water depths were overestimated across all models with RMSE between 0.63 and 0.84 m and MAE between 0.51 and 0.63 m (Figure 8). The SVR model outperformed the other models with RMSE = 0.63 m and MAE = 0.51 m. Overestimation across models is likely attributed to the differences between the measured data used to train models and the test data. This is because the trained models have a higher tendency to predict depths in the range of the training data which reduced their performance on data with different characteristics (

Table 7. Descriptive statistics of the training, testing and transfer data to show how characteristics of the testing and transfer sites compare to the training sites (values in m).

Statistic	Training Data	Manual Data	HOBO Test 1	HOBO Test 2	P36	G-3272	NP201
Mean	0.73	0.12	0.21	0.50	0.34	0.18	0.30
Median	0.55	0.00	0.01	0.22	0.28	0.13	0.23
Std. Dev	0.81	0.17	0.35	0.81	0.21	0.18	0.28
Min	0.00	0.00	0.00	0.01	0.04	0.00	0.00
Max	5.26	0.61	1.42	3.23	0.76	0.56	0.88
25th percentile	0.25	0.00	0.01	0.02	0.18	0.00	0.08
75th percentile	0.93	0.24	0.36	0.49	0.48	0.38	0.45

).

Further testing was performed utilizing the data from the two HOBO loggers to corroborate the assertion that the target test site should have a comparable hydrologic regime. HOBO test site 2 had more similar characteristics to the training data than test site 1 (Table 7). Figure 9a shows that all the models overestimated the low to dry conditions but could capture peak depths well at the HOBO test site 1 despite differences in hydrologic regimes compared to training data. A similar performance was seen at the HOBO test site 2 (Figure 9b). Models performed better on the HOBO logger data which had similar characteristics to training data (

Table 8. Performance of SDB models at testing locations where measurements were taken manually and at two HOBO logger locations in West Tennessee.

Site	Model	RMSE (m)	MAE (m)	${\mathit{R}}^2$
Manually-measured sites	Stumpf	0.79	0.63	0.55
	Bramante	0.83	0.62	0.58
	RF	0.92	0.64	0.43
	SVR	0.63	0.51	0.55
	KNN	0.74	0.59	0.55
HOBO Test Site 1	Stumpf	0.27	0.24	0.65
	Bramante	0.22	0.20	0.67
	RF	0.28	0.26	0.49
	SVR	0.25	0.24	0.65
	KNN	0.26	0.25	0.70
HOBO Test Site 2	Stumpf	0.42	0.29	0.88
	Bramante	0.35	0.27	0.87
	RF	0.47	0.32	0.79
	SVR	0.51	0.31	0.83
	KNN	0.40	0.28	0.91

). The results suggest models may be transferable to sites with similar environmental characteristics and hydrologic regimes. We suggest further evaluation at similar sites in the Mississippi alluvial plains to test this hypothesis. Despite superior performance at the HOBO test sites, no model was able to predict dry conditions because training data did not include water depths of zero (Figure 9).

3.4. Model Transferability

Results obtained at transfer sites resembled those from test sites. All models generally overestimated water depth at transfer sites (Figure 10). However, despite slight overestimation, the models successfully captured temporal water depth fluctuations across all the sites, aligning with the results obtained at the test sites. Notably, superior model performance was evident at transfer site P36, with RMSE ranging from 0.15 to 0.29 m and MAE from 0.12 to 0.60 m across all models, in comparison to transfer sites G-3272 and NP201 (

Table 9. Performance of SDB models at transfer locations: P36, G-3272, and NP201 in the Everglades Depth Estimation Network (EDEN), Florida.

Site	Model	RMSE (m)	MAE (m)	${\mathit{R}}^2$
P36	Stumpf	0.19	0.18	0.16
	Bramante	0.15	0.58	0.11
	RF	0.29	0.12	0.23
	SVR	0.16	0.43	0.12
	KNN	0.19	0.46	0.15
G-3272	Stumpf	0.50	0.41	0.46
	Bramante	0.46	0.85	0.42
	RF	0.48	0.30	0.41
	SVR	0.35	0.62	0.31
	KNN	0.43	0.70	0.40
NP201	Stumpf	0.91	0.73	0.86
	Bramante	0.94	0.67	0.83
	RF	0.88	0.45	0.73
	SVR	0.61	0.61	0.54
	KNN	0.82	0.66	0.69

). The distribution of the residual errors at P36 generally fell within −0.2 to 0.2 m, with few outliers. The Bramante and SVR models outperformed the other models, as illustrated in Figure 11a. This strong performance at P36 can be attributed to the similarity between the observed data and the training data, as indicated in Table 7.

The characteristics of the observed data at NP201 were similar to the training data. However, the performance of the models at site G-3272 was better than what was obtained at NP201. This observation was contrary to the trends seen during testing, where the models performed better at sites with observed data having similar characteristics to the training data. A review of the site characteristics of NP201 indicated the presence of sawgrass and emergent marsh vegetation, in contrast to the open water sites from which the training data were collected. The presence of this vegetation affects the spectral reflectance from these sites and may account for the poor performance of the models. Similar performance was seen at the manually measured test sites where the data were collected from flooded corn fields. The results obtained at the transfer sites demonstrated the potential transferability of the models. The model parameters may also be useful for estimating the parameters for the same models in different areas.

The test results indicate that the methods employed are constrained by environmental and hydrologic characteristics. The methods cannot be applied to forested or heavily vegetated wetlands which prevent satellites from detecting the reflectance from the bottom of surface water bodies. Longer satellite revisit times and cloud contamination may affect the ability to capture important hydrologic events.

4. Conclusions

This work presents empirical and machine-learning methods for estimating water depths of inland wetlands. The Lyzenga [32] and Stumpf et al. [34] empirical bathymetric models were modified to extend their application from clear coastal waters to turbid inland wetlands, which outperformed machine learning models (i.e., Random Forest, Support Vector Regression, and k-Nearest Neighbor) based on selected performance evaluation metrics. All the models developed performed poorly in deeper waters with depths exceeding 3.0 m. The inclusion of the level of turbidity of the waters improved the performance of the machine learning models in deep waters but not the empirical models. The effect of turbidity in shallow waters was found to be insignificant because estimates in shallow waters are largely controlled by bottom reflectances and not water column properties. Models demonstrated promising transferability, suggesting they may be transferred and generalized to other locations with similar environmental and hydrologic characteristics; however, we suggest additional testing and research are required. Applying these models to a satellite image of a wetland allows practitioners the ability to estimate whole-wetland scale hydrology from one or a few measured locations.

Despite predictable limitations, this work is a significant step forward in accurately estimating water depth in inland wetlands using remote sensing. Therefore, its contribution and applications for wetland monitoring, management, and conservation planning are significant and numerous. For example, inland wetlands are vast and store 10-fold more carbon than coastal wetlands [9]; therefore, natural resource economists may consider using remotely-sensed bathymetry models to accurately assign economic wetland valuations (i.e., carbon sequestration potential, stormwater retention, and wildlife values; [3,17]). Depth is challenging to incorporate into conservation planning, despite known water levels that promote multi-species abundance and diversity [11,52]. Indeed, wetland depth was excluded from the Lower Mississippi Alluvial Valley Migratory Bird Joint Venture’s (JV) conservation planning models because water depth information was insufficient and monitoring was impractical (North American Waterfowl Management Plan [NAWMP] [11,53,54]). Therefore, conservation planners should consider using remotely-sensed bathymetry models to monitor wetland conditions across broad geographic extents needed for effective management of migratory waterbirds (NAWMP [53]; 2018; [55,56,57]). Many other applications exist across wetland conservation, management, and policy domains. Undoubtedly, the presented remote sensing and depth estimation models will empower practitioners to make more informed decisions (e.g., prioritize protection, and optimize management), at reduced cost, and across larger geographic extents.

5. Recommendations

Future work should continue to test the transferability of remotely-sensed bathymetry models in regions outside of the Mississippi alluvial plains and evaluate these and other integrated algorithms to estimate depth in forested wetlands. Lastly, refinement of machine-learning models to estimate water depth is prudent; these methods appeared to perform well in this study in deeper wetlands where empirical algorithms struggled.

Rapid development and research using remote sensing technologies for wetland monitoring is encouraging but limitations remain in this application and elsewhere (e.g., [30,58]). For example, remotely-sensed bathymetry is difficult to accurately quantify in forested wetlands or other densely vegetated areas (i.e., coastal marsh). Cloud contamination and longer satellite revisit times of the Landsat imagery used in this study can also exacerbate the problems in using remotely sensed products in wetland monitoring. It is recommended that future works consider using Sentinel-1 imagery to solve the challenges due to cloud contamination and satellite revisit times and also explore other depth estimation techniques that can be applied to forested wetlands.