Deriving Tidal Structure From Satellite Image Time Series

In shallow coastal regions, tides often control the water flux, which in turn directs sediment transport, nutrient delivery, and geochemical gradients. However, tides in shallow areas are spatially heterogeneous, making it challenging to constrain the geographic structure of tidal phase and amplitude without extensive networks of tide gauges. We present a simple remote sensing method for deriving tidal structure from satellite time series. Our method is based on two observations: (1) Tidally driven variations in water depth can be detected as changes in pixel intensity in optical satellite imagery, and (2) repeating passes by an orbiting satellite capture a region at different phases of the tidal cycle. By stacking multiple satellite acquisitions of a shallow bank, we can compute the relative tidal phase and amplitude for each pixel location, thereby resolving a detailed map of tidal propagation and attenuation. While our method requires a set of local water‐depth measurements to calibrate the color‐to‐depth relationship and compute tidal amplitude (in meters), our method can produce spatial estimates of tidal phase and relative amplitude without any site‐specific calibration data. As an illustration of the method, we use Landsat imagery to derive the spatial structure of tides on the Great Bahama Bank, estimating tidal phase and amplitude with mean absolute errors of 15 min and 0.15 m, respectively.


Introduction
The behavior of oceanic tides on shallow continental shelves-and the implications for topics as broad as the evolution of the Earth-Moon orbit (Coughenour et al., 2009;Hansen, 1982) and the global mixing of the ocean (Munk & Wunsch, 1998)-has occupied the minds of scientists for over a century (Jeffreys & Shaw, 1920;Taylor & Shaw, 1919). The advent of satellite altimetry in the early 1990s revolutionized our understanding of global oceanic tides (Cartwright & Ray, 1990;Egbert & Ray, 2000;Tierney et al., 2000). However, the coarse track spacing (∼100 km) of satellite altimeters means that we depend on models-which assimilate satellite altimetry and tide gauge observations-to obtain continuous tidal predictions. While global tide models can constrain tides within 2 cm in the open ocean (Shum et al., 1997), bathymetric uncertainty and coarse grid spacing cause these models to break down in shallow regions such as the Great Bahama Bank (GBB), where tidal lags can exceed 4 hr over distances as small as 20 km ( Figure 1).
Three-dimensional numerical circulation models (e.g., Fringer et al., 2006;Hervouet, 2000;Luo et al., 2013), which now are applied routinely to embayments (Weisberg & Zheng, 2006), banktops (Purkis et al., 2017(Purkis et al., , 2019, and estuaries (Zheng & Weisberg, 2004), are a powerful tool for understanding water circulation in shallow areas. However, these models typically require bathymetric, wind, and tide data for model initialization and forcing. Available remote sensing and weather reanalysis products may be used to constrain some aspects of these hydrodynamic models. For example, satellite scatterometer data constrain surface winds, and water depth can be derived from satellite imagery in shallow waters with low turbidity (e.g., Geyman & Maloof, 2019;Lyzenga, 1978Lyzenga, , 1981Stumpf et al., 2003). Unfortunately, the geographic structure of tidal phase and amplitude often is unknown a priori, and the paucity of tide gauge data often forces the use of low-resolution global tide models in hydrodynamic simulations (e.g., Purkis et al., 2017Purkis et al., , 2019. We present a simple method for extracting tidal structure from satellite time series. Our method is based on two observations. First, tidally driven variations in water depth can be detected as changes in pixel intensity  (Broecker & Takahashi, 1966), tidal lags can exceed 4 hr over distances of ≤20 km. (a) A map of the region NW of Andros Island with the tidal lags (in hours) at five tide gauges reported relative to the NOAA tide predictions station at Fresh Creek, Andros Island. "TOTO" stands for "Tongue of the Ocean," the ∼1to 2-km-deep trough east of Andros Island. The gray lines represent the positions of ∼300,000 echosounder depth measurements (Geyman & Maloof, 2019a). Only the ∼50,000 depth measurements acquired within ≤5 km of a tide gauge, where our understanding of the tides is well constrained, are used as calibration data. in optical satellite imagery. Second, repeating passes by an orbiting satellite capture a region at different phases of the tidal cycle. Thus, by stacking multiple satellite acquisitions of a shallow bank, we can compute the relative tidal phase and amplitude for each pixel location, thereby resolving a detailed map of tidal propagation and attenuation.

RESEARCH ARTICLE
As an illustration of the method, we use Landsat imagery to derive the geographic structure of tides on a portion of the GBB. We choose the GBB as a proof-of-concept study because shallow water carbonate platforms represent the most important ancient archive of sea level, environmental, and ocean geochemical information, so studying modern analog settings such as the Bahamas is a key component of any translation of the rock record into paleoenvironmental interpretations (e.g., Broecker & Takahashi, 1966;Dyer et al., 2018;Geyman & Maloof, 2019b;Harris et al., 2015;Maloof & Grotzinger, 2012;Oehlert et al., 2012;Patterson & Walter, 1994;Swart et al., 2009). It has become clear in modern settings that tidal pumping of water and nutrients plays a critical role controlling the transport and chemistry of carbonate sediment, but global tide models still cannot resolve tidal dynamics on shallow platforms (Shum et al., 1997).

Study Area and Field Observations
The shallow waters NW of Andros Island on the GBB are well suited for study of spatial variability in tidal structure, since tidal amplitude can vary by a factor of 2 and tidal phase can shift by ∼4 hr over distances of ≤20 km ( Figure 1). Banktop tides are predominantly semidiurnal, with strong spring-neap amplitude fluctuations (Figures 1c and 4a). We use a network of five tide gauges, each recording 5 to 18 days of data in June 2017 and June 2018 (Figures 1d-1h), as ground-truth to evaluate the results of our model. However, only one reference tide gauge is required to implement our method. We use the NOAA tide predictions station at Fresh Creek, Andros Island ( Figure 1) as our reference. Note that it is not necessary for the reference tide gauge to be located within the extent of the imagery. The singular requirement is that the tidal character (i.e., the relative mixture of tidal constituents; Figure 1c) at the reference tide gauge is similar to that in the region of interest.
To calibrate the conversion from satellite image color to water depth (see section 2.3), we use ∼50,000 depth measurements collected with a BioSonics MX Aquatic Habitat Echosounder (single frequency, 204.8 kHz) (Geyman & Maloof, 2019a).

Inferring Tidal Structure From Satellite Time Series
Landsat images of the GBB always are acquired at approximately the same time-between 15:20 and 15:40 GMT. However, because tides on the GBB are dominated by the lunar semidiurnal (M2) tide (Figure 1c), the high tide is pushed backward by about 50 min each day. Thus, over time, Landsat images capture the shallow bank at all phases of the tidal cycle ( Figure 2).

Estimating Water Depth From Multispectral Satellite Imagery
Light is attenuated exponentially as it travels through the water column, causing deep water to appear darker than shallow water. Therefore, the color variation in a stack of satellite images records the propagation of tides across a shallow bank. To compute changes in tidal amplitude (in meters) across the banktop, we must convert pixel intensity to water depth. Lyzenga (1978Lyzenga ( , 1981 proposed a basic algorithm for calculating water depth from satellite imagery, and applications and adaptations of the original Lyzenga algorithm have proven widely successful in shallow waters around the world (Bierwirth & Burne, 1993;Dierssen et al., 2003;Kerr & Purkis, 2018;Lyzenga et al., 2006;Mumby et al., 1998;Philpot, 1989;Polcyn et al., 1970;Stumpf et al., 2003).
We use a simple modification of the Lyzenga (1978) linear algorithm in which we segment the satellite image into zones of spectral homogeneity and calibrate different color-to-depth relationships for each zone before applying a fuzzy k-nearest neighbor fit to solve for the final bathymetry. This cluster-based regression (CBR) method of Geyman and Maloof (2019a) reduces error by ≥40% relative to previous bathymetric estimates for the region in Figure 1 (Geyman & Maloof, 2019a), since it better accounts for variability in the spectral properties of the sea floor. When applied to the 56 Landsat images in our compilation, the CBR method produces depth estimates with mean absolute errors (MAEs) of 0.14-0.26 m and R 2 values of 0.70-0.88 ( Figure 3; Table S1).

Correcting for Spring-Neap Amplitude Fluctuations
The 56 Landsat images in our compilation sample all portions of the tidal cycle ( Figure 2). Of course, each Landsat acquisition captures a different tidal cycle, since the revisit period of Landsat satellites is 16 days. Tidal amplitude on the GBB varies by nearly a factor of 2 over the course of a typical bimonthly spring-neap cycle ( Figure 4a). Therefore, to facilitate comparison between images captured from neap to spring tide, we scale all Landsat-derived estimates of tidal amplitude byĀ FC ∕A FC , whereĀ FC is the mean tidal amplitude at Fresh Creek (0.93 m) and A FC is the contemporaneous tidal amplitude at Fresh Creek. For example, if a particular Landsat-derived tidal estimate is +0.3 m for a given location on the bank, but the tide station at Fresh Creek indicates that the tidal amplitude on the day of image acquisition was just 0.50 m, the tidal offset observation at the pixel of interest will be scaled by 0.93/0.50.   (Table 1) reflect the comparison between tide gauge records and ∼200 random pixel locations within 500 m of each tide gauge. To prevent biasing our error estimates, we exclude bathymetric data acquired near a given tide gauge when estimating error for that particular tide gauge. For example, in (a), which represents model predictions for the Joulters (2017) tide record, only bathymetric data acquired near the other four tide gauges were used to calibrate the CBR algorithm (section 2.3). The gray "+" marks, and "x" marks indicate depth estimates from Landsat 5 and Landsat 8 acquisitions, respectively. The text box shows the phase and amplitude error for each sinusoidal fit (colored line).

Computing Tidal Phase and Amplitude
In regional and global tide models, tides typically are defined in terms of the amplitudes (meters) and phases (radians or degrees) of a set of tidal constituents, such as the lunar semidiurnal (M2) and solar diurnal (S1) tides (Lyard et al., 2006). Our 56 sets of Landsat-derived water depth estimates are too sparse (and temporally incoherent) to compute the phase and amplitude parameters of the prominent tidal constituents for each pixel. Therefore, we describe GBB tides in terms of a phase lag and amplitude change from the tides at Fresh Creek ( Figure 4).  Figure 2 summarizes our method of inferring tidal phase and amplitude either (A) with or (B) without local water depth data. If local water depth data are available, we use the CBR algorithm to generate a bathymetric map for each of the 56 Landsat images (section 2.3). Thus, for every pixel location on the bank, we obtain a vector of 56 water depth estimates (and their associated uncertainties) and a corresponding vector of positions in the tidal cycle (measured as hours from high tide at Fresh Creek). We use the Nelder-Mead multidimensional unconstrained nonlinear minimization algorithm (Nelder & Mead, 1965) to solve for the best fit sine wave with a 12.42-hr period. Note that 12.42 hours is not only the period of the dominant tidal constituent at Fresh Creek (M2, Figure 1c) but also the mean tidal period observed in 13 years of high/low tide data for Fresh Creek ( Figure S6). For each pixel location, the sine wave's phase offset from the Fresh Creek reference provides an estimate of tidal lag, and the sine wave's amplitude records the local tidal amplitude. If local bathymetric data are not available, we cannot use the CBR algorithm to map spectral reflectance to water depth. However, our method still can be used to generate estimates of tidal phase and relative amplitude (Figure 2).

Model Validation
To evaluate the accuracy of our tide model, we compare the predictions of tidal phase and amplitude to the measured values at each of our five gauges (Figure 1). Note that we use the tide gauge records to correct bathymetric data used in the CBR algorithm (section 2.3). Therefore, to ensure a fully independent evaluation of model accuracy, we iteratively train the tide model using bathymetric data adjacent to four of the five tide gauges. We then compare the estimated phase and amplitude at the location of the held-out tide gauge to the ground-truth measurements. Figure 5 compares the model-derived tide estimates at five nearby pixel locations to the true tidal phase and amplitude. The tide model that incorporates local bathymetric data Figure 7. Sensitivity of the estimated tidal phase and amplitude to k, the number of spectral clusters used in the CBR algorithm (Geyman & Maloof, 2019a). For k ≥ 2, all values of k yield similar predictions for tidal phase (a-f) and amplitude (g-l).
yields MAEs for tidal phase and amplitude of 0.25 hr and 0.15 m, respectively. The tide model that does not rely on any local ground-truth data produces a larger average error of 0.78 hr ( Table 1).
The typical scatter of ≤0.3 m about the sinusoidal fits in Figure 5 is unsurprising given that the bathymetric model, which maps image color to water depth (Figure 3), has a MAE of 0.14-0.26 m (Table S1). However, some of the noise in the records of tidal position versus water depth may be caused by real changes in bathymetry during the 14-year Landsat time series. The high-energy Joulter Cays region, which experiences water depth changes on diurnal to annual timescales caused by migrating subaqueous dunes and ooid shoals (Figure 6), produces noisier records than the other sites ( Figure 5). The 30-m resolution of the Landsat images is too coarse to resolve small sand bodies (Figure 6b), but time series analysis of the residuals between the individual depth observations and the sinusoidal fits in Figure 5 shows patterns consistent with migrating ooid shoals (Figures 6c-6e). Thus, the residuals in tidal amplitude may provide a window into understanding the timescales of bedform migration and sediment flux. Moreover, while the elevated noise in the Joulters sites causes the tidal amplitude to be underestimated, the tidal lags at those sites are still correct to within ≤0.19 hr (11 min) (Table 1).

Model Sensitivity
The only tunable parameter in our tide model is k, the number of spectral clusters used during the conversion of reflectance to water depth in the CBR algorithm. Beyond k = 1, all models converge on similar tidal phase and amplitude estimates (Figure 7). Therefore, because Geyman and Maloof (2019a) showed that models with lower k values perform better when extrapolating from limited calibration data sets, we use k = 2. Figures 8a and 8b show the tide model's phase and amplitude predictions for the GBB. The lag relative to Fresh Creek is close to zero along the margin closest to the Tongue of the Ocean (Figure 1). The shallow waters on the west side of Andros Island and Joulter Cays permit long lags (∼4 hr) and small tidal amplitudes (∼0.3 m). If no local bathymetric data were available, we could generate similar predictions of tidal lag and relative amplitude (Figures 8c and 8d). The two predictions of tidal lag typically differ by ≤1 hr ( Figure S4).

Uncertainty Analysis
To understand the uncertainty associated with the tide estimates in Figure 8, we apply a bootstrap resampling approach in which, for each pixel, we fit a sine wave to a random sample of 70% the data (40 of the 56 Landsat-derived depth estimates). We repeat this procedure with 20 random subsamples and then report the interquartile range (25th to 75th percentile) in tidal lag and amplitude. Figure 9 illustrates that shallow Figure 8. Spatial structure of tidal phase (a) and amplitude (b) generated by the tide model that uses local water depth data ( Figure 2A). Notice that the majority of the phase lag and amplitude attenuation occurs within ∼10 km of the shelf break where the GBB gives way to the Tongue of the Ocean (Figure 1). If no local bathymetric data were available, we could generate similar predictions of lag (c) and relative amplitude (d). The phase estimates in (a) and (c) typically differ by ≤1 hr (Supplementary Information, Figure S4).
waters support the smallest uncertainties in both tidal lag and amplitude, which likely is because the CBR algorithm performs best in shallow waters (Geyman & Maloof, 2019a).

Discussion
Regions such as the GBB, with very shallow waters and significant restriction from the global ocean, have eluded previous attempts at modeling tides from remote sensing data. On the GBB, the gradients in tidal phase are too steep to resolve directly from satellite altimetry missions, which often have track spacing ≥100 km (∼315-, ∼140-, and ∼80-km spacing at the equator for TOPEX/Poseidon, JASON-1, and ERS5 1, respectively). Tide models assimilating satellite altimetry data with the barotropic shallow-water equations provide tidal estimates with complete banktop coverage. However, coarse grid spacing and large relative bathymetric errors cause these models to be inaccurate in shallow waters (Shum et al., 1997). For example, the high-resolution (0.125 • × 0.125 • ) global tide model FES2014 (Carrere et al., 2015) predicts 0 min of tide lag and no amplitude attenuation between Joulter Cays and Triple Goose Creek (TGC), where we document 4.5 hours of lag and 45% reduction in tidal amplitude (Figure 4). We quantify uncertainty as the interquartile range of lag (a) and amplitude (b) estimates obtained by iteratively applying the tide model to random subsets of the Landsat data set (section 3.3). Both tidal lag and amplitude are poorly constrained for the deeper waters along the western margin of the study area, which may reflect the fact that the conversion from pixel intensity to water depth is less accurate in deeper waters ( Figure S5). (c) Bathymetry from Geyman and Maloof (2019a).
The fact that global tide models, which now boast errors of <2 cm for the open ocean (Shum et al., 1997), remain inaccurate for shallow waters is important for two reasons. First, accurate representation of banktop tides is necessary for forcing hydrodynamic models, which enable understanding of the link between water flux, geochemical gradients, sediment transport, and biosedimentary facies. Even our tidal estimates that require no site-specific calibration data (Figures 8c and 8d) are capable of representing the banktop tidal propagation and attenuation missed by global models such as FES2014. Second, shallow shelves are significant contributors to the global budget of tidal dissipation (Egbert et al., 2004). Our simple method for estimating tidal structure could be scaled up and applied to coarser optical imagery such as MODIS in order to quantify tidal dissipation in shallow water regions across the globe.
Although freely available Landsat images yield empirical tidal estimates with continuous coverage and 30-m resolution, our optical method suffers from a number of disadvantages compared to the altimetric method. First, optically derived tidal estimates are sensitive to changes in water turbidity, which may vary with tidal phase, potentially biasing the retrieved tidal signal and restricting application to clear (nonturbid) waters like those of the GBB. Second, optical tidal estimates are susceptible to changes in local bathymetry and bottom reflectance. Real changes in benthic character over the 14-year Landsat time series (Figure 2) may violate the assumption that changes in each pixel's reflectance represent tidally driven changes in water depth ( Figure 5) and may contribute to the uncertainty of our final estimates. However, new laser altimetry data from ICESat-2 (Abdalati et al., 2010), which can provide elevation estimates for both the sea surface and sea bottom (Forfinski-Sarkozi & Parrish, 2016), could serve as a cross-calibration for optically derived tidal estimates. Moreover, the optical method's sensitivity to seabed properties also represents one of its greatest strengths, since our approach can be used to quantify sediment transport ( Figure 6) and changes to benthic ecology on decadal time scales.