Mapping Water Surface Elevation and Slope in the Mississippi River Delta Using the AirSWOT Ka-Band Interferometric Synthetic Aperture Radar

AirSWOT is an airborne Ka-band synthetic aperture radar, capable of mapping water surface elevation (WSE) and water surface slope (WSS) using single-pass interferometry. AirSWOT was designed as a calibration and validation instrument for the forthcoming Surface Water and Ocean Topography (SWOT) mission, an international spaceborne synthetic aperture radar mission planned for launch in 2022 which will enable global mapping of WSE and WSS. As an airborne instrument, capable of quickly repeating overflights, AirSWOT enables measurement of high frequency and fine scale hydrological processes encountered in coastal regions. In this paper, we use data collected by AirSWOT in the Mississippi River Delta and surrounding wetlands of coastal Louisiana, USA, to investigate the capabilities of Ka-band interferometry for mapping WSE and WSS in coastal marsh environments. We introduce a data-driven method to estimate the time-varying interferometric phase drift resulting from radar hardware response to environmental conditions. A system of linear equations based on AirSWOT measurements is solved for elevation bias and time-varying phase calibration parameters using weighted least squares. We observed AirSWOT WSE uncertainty of 12 cm RMS compared to in situ water level measurements when averaged over an area of 0.5 km2 at incidence angles below 15◦. At higher incidence angles, the observed AirSWOT elevation bias is possibly due to residual phase calibration errors or radar backscatter from vegetation. Elevation profiles along the Wax Lake Outlet river channel indicate AirSWOT can measure WSS over a 24 km distance with uncertainty below 0.3 cm/km, 8% of the true water surface slope as measured by in situ data. The data analysis and results presented in this paper demonstrate the potential of AirSWOT to measure water surface elevation and slope within highly dynamic and spatially complex coastal environments.


Introduction
Quantifying and monitoring the volume and distribution of water is important, particularly as climate change and other factors cause additional stress to Earth's limited water resources [1,2]. Accurate measurements of these water resources are important in planning for both overabundance and dearth of water, as in floods and droughts, respectively. Rivers and lakes are particularly important due to their ease of access and role in flood hazards. Unfortunately, in situ data regarding surface water //deltax.jpl.nasa.gov/) will use several airborne remote sensing instruments, including AirSWOT, in order to study the relative contribution of vegetation and channel-network to soil accretion rates. We also seek to demonstrate AirSWOT's capabilities to generate accurate WSE and WSS measurements supporting potential implementations of the Delta-X framework to other study areas.
This paper is structured as follows. We introduce the data used in this study in Section 2 and present an overview of the AirSWOT instrument in Section 2.1. The in situ data collection and processing is discussed in Section 2.2. We introduce the operational interferometric phase calibration methodology in Section 3.1, followed by the proposed calibration approach adapted to these coastal landscapes. In Section 3.2 we describe the filtering and averaging methodology we used to estimate WSE from the noisy AirSWOT data. The results are presented in Section 4. We present and discuss the AirSWOT WSE results in Section 4.1. In Section 4.2, we show relative WSE estimation over time. AirSWOT-derived profiles of WSE and WSS in the Wax Lake Outlet river channel are shown in Section 4.3. Our conclusions are discussed in Section 5.

AirSWOT Overview
AirSWOT consists of three main instruments. The primary instrument is called the Ka-band SWOT Phenomenology Airborne Radar, or KaSPAR. KaSPAR is the interferometric synthetic aperture radar used by AirSWOT to measure surface water levels. The other two instruments are an Applanix POSAV 610 integrated GNSS and IMU (global navigation satellite system and inertial measurement unit), and a Cirrus color infrared digital camera system. The GNSS and IMU are used to precisely track the aircraft position and remove the effects of aircraft motion and attitude on the radar data. The digital camera can be used to map the extents of surface water in cloud-free conditions (e.g., [37]). In this study the camera was not installed, and instead, we used a pre-existing land cover mask generated from Sentinel-2 data by Thomas et al. [38] to delineate land and water.
KaSPAR uses six different antennas to collect interferometric data in a variety of imaging modes and with different cross-track and along-track spatial baselines. In this study, we focus on data from the KaSPAR outer swath mode, which has incidence angles ranging nominally between 4-25 • . In practice, large errors are often observed at the edges of this range, such that the usable swath is generally smaller (see Section 4.1).
The main characteristics of AirSWOT and the AirSWOT outer swath data are shown in Table 1. In addition to the outer swath mode, KaSPAR has an inner swath mode with incidence angles ranging between 0-7 • , designed to more closely reflect the range of incidence angles that will be collected by SWOT. Compared to the outer swath mode, the inner swath mode uses a smaller cross-track antenna beamwidth and a higher range bandwidth of 450 MHz, yielding a range resolution of ≈ 40 cm (compared to ≈ 2 m for the outer swath). However, for this study, we focus exclusively on the outer swath mode which, advantageously, covers a wider surface swath.
AirSWOT uses both cross-track and along-track interferometry to accurately measure water surface elevation. Cross-track interferometry is sensitive to the vertical position of radar scatterers in a scene [14], allowing cross-track interferometry to be used to estimate elevation. However, using cross-track interferometry alone to estimate water surface elevation would result in biased WSE estimates. This is because water is a moving target, and the velocity of water in the direction of the radar line of sight will cause a doppler shift, which changes the interferometric phase [39,40]. Along-track interferometry can be used to estimate the motion of the water scatterers [39,40]. AirSWOT uses along-track interferometry to estimate the doppler shift introduced by water motion, and then compensate for the effects of that doppler shift in the cross-track interferograms. The cross-track and along-track interferometric data are combined into a motion-compensated cross-track interferogram, from which the water surface elevation can be estimated with the effects of the motion-induced interferometric phase bias removed. Since KaSPAR is a single-pass multi-static InSAR system, the cross-track and along-track data are collected simultaneously. For the outer swath mode, KaSPAR transmits using one antenna and receives using three antennas (one of which is the same as the transmit antenna), forming two interferometric baselines. The first baseline is in the cross-track direction, while the second baseline is primarily along-track, with a small cross-track component. More details on the KaSPAR instrument can be found in [23].

In Situ Water Level Data
We used in situ water level measurements to calibrate and validate the AirSWOT estimates. Several stations within the study regions are managed by the United States Geological Survey (USGS), the National Oceanic and Atmospheric Administration (NOAA), and the Coastwide Reference Monitoring System (CRMS) [41]. These organizations publicly distribute their data online. We also used a Louisiana State University (LSU) dataset of in situ water level collected within the Wax Lake Delta, on Mike Island. A summary of the USGS, NOAA, CRMS, and LSU water level stations used in this study are shown in Table 2. In the table, the first two stations (Calumet and Caillou Bay) were used in the AirSWOT calibration procedure (Section 3.1). The remaining stations were used as independent validation of the AirSWOT water level estimates. In the table, the approximate geographic position of each station is listed. For this study, we used NAVD88 as the common vertical datum to compare between the remote sensing and in situ data. In the last column of the table, we list how the in situ water levels with respect to the gauge datums were converted to the NAVD88 vertical datum (using the GEOID12B geoid model). The USGS stations were easy to convert, as each USGS station used in this study had a published offset with respect to NAVD88. There were two USGS stations within our study area (at Crewboat Channel and Caillou Lake) which did not have publicly available datum information (or, in the case of Caillou Lake, had datum information using an unknown geoid model). These stations were excluded from the analysis. While the Calumet station has a published vertical offset of −17.07 cm with respect to NAVD88, this did not explain a bias we initially observed between AirSWOT and the station. After inquiring with the USGS Baton Rouge Office we learned the datum offset was with respect to an older geoid model (GEOID99) when the Calumet station was last surveyed in 1996 [42]. We converted the Calumet water levels from GEOID99 to GEOID12B using the publicly available NOAA VDatum software tool [43] and determined an additional offset of −20.7 cm.
Water levels from the two NOAA water level stations are published with respect to local tidal datums and were converted to NAVD88 using NOAA VDatum [43]. VDatum reports the uncertainty associated with this datum conversion, which was 15.7 cm for both the Eugene Island and Amerada Pass stations. The LSU Mike Island water level gauge had its datum surveyed using real-time kinematic (RTK) positioning in February of 2017. We assume sinking of the water level gauge to be negligible between the gauge and AirSWOT survey dates.
The CRMS stations report their water levels in NAVD88, with the geoid model used for each water level sample included in the data. All of the CRMS data used in this study were already in the GEOID12B geoid model, therefore no vertical datum conversion was required for these data.
A first assessment of AirSWOT WSE estimates revealed conflicting behaviours between CRMS and non-CRMS stations. While the USGS and NOAA stations tend to be located in large river channels or in large areas of open water, many CRMS stations are located in small distributary channels often located near emergent wetland vegetation. The LSU station was located near the middle of Mike Island, which was flooded during the time period of this study. Therefore, we consider separately the non-CRMS and CRMS stations throughout this paper, particularly when discussing the results in Section 4. Figure 1 shows an overview map of the study area, with the locations of the in situ water level stations superimposed on a map of the calibrated AirSWOT elevation data, which should not be confused with a map of WSE. The map represents the elevation estimated from the AirSWOT interferometric phase, which is biased compared to the true WSE in vegetated areas. In order to estimate WSE from the AirSWOT elevation products, additional filtering and averaging is required (Section 3.2). The stations are colored according to their source. The two stations used as GCPs for AirSWOT calibration (Calumet and Caillou Bay) are annotated with their names. We chose to use Calumet and Caillou Bay as the calibration stations because they are located near the corners of the study area and are representative of the high and low ends of the observed water level range, respectively. Figure 2 shows a map zoomed in more closely on the Wax Lake and Atchafalaya River deltas. In this map, some of the stations used for validation (Avoca Island Cutoff, Mike Island, Eugene Island, Amerada Pass, and CRMS 0479) are annotated with their names.
Unfortunately, the NOAA stations at Eugene Island and Amerada Pass were located outside the bounds of the AirSWOT data coverage. For the purposes of our analysis, we treated these stations as if they were located within the coverage area of nearby AirSWOT flight lines by creating a "virtual location" for each station. In Figure 2, the true location of each NOAA station is shown in dark blue, while the virtual location used in our analysis is shown in light blue. For Amerada Pass, the virtual location is approximately 2.3 km from the true position, and has been moved from a small deltaic river channel into the open water outside the delta. This is also helpful because vegetation near river channels can introduce biases into the AirSWOT water level estimates (as discussed in Section 4.1). The distance between the true and virtual locations for the Eugene Island station is 3.2 km, with the virtual location lying inside the nearest AirSWOT swath. We justify this by noting that the Eugene Island and Amerada Pass stations are located in the Atchafalaya Delta, close to the ocean, where it is expected the WSS is relatively flat compared to the steeper upstream WSS. To determine the potential WSE error introduced by moving these stations, we calculated the WSS between the Eugene and Amerada stations from the in situ data, finding that the WSS varied between 0.04 and 0.88 cm/km during the AirSWOT acquisition period spanning 3 days. By moving the Eugene Island station by 3.2 km, the WSS values may result in a change in WSE between 0.13 and 2.83 cm. For the Amerada Pass station, the change in WSE is between 0.09 and 2.02 cm. This change in WSE between the original and shifted station locations is significantly smaller than the elevation uncertainty of the stations' vertical datum conversion (15.7 cm) or the AirSWOT data itself (as shown in Section 4).

AirSWOT Phase Calibration Procedure for Coastal Wetlands
Phase calibration is a critical part of InSAR processing [44][45][46][47]. In the case of AirSWOT, the phase calibration consists of estimating unintended differential phase delays between radar channels in the KaSPAR instrument. These differential phase delays cause an InSAR phase bias if not removed. An InSAR phase bias introduces a cross-track slope in the measured elevations, because the relationship between interferometric phase and elevation varies with incidence angle [14]. For AirSWOT, this phase bias varies with time-due in part to changes in instrument temperature-introducing a time-varying cross-track elevation slope if not estimated and removed. This phase drift was found to differ between flights and must be estimated from the InSAR data itself.
The KaSPAR instrument contains internal calibration loops that were designed to measure the phase delays between certain points within the radar hardware. Unfortunately, in practice these calibration loops are insufficient to fully characterize the observed phase drift, leaving residual phase calibration errors. The AirSWOT processor therefore uses the calibration loops to provide an initial estimate of the differential phase delays, which are then refined using data-driven methods. Note that due to hardware and platform differences between KaSPAR and KaRIn (Ka-band Radar Interferometer, the primary instrument on SWOT), the phase calibration methods we discuss here are only applicable to KaSPAR. For information on SWOT calibration plans, see [48].
The operational AirSWOT phase calibration procedure fits a set of phase correction coefficients for each AirSWOT flight line: a constant phase bias and, optionally, a phase bias term that is a linear function of along-track distance. The linear phase bias term is only used for flight lines where the phase correction changes significantly over the course of the flight line. The phase corrections usually change rapidly at the beginning of each flight, then gradually become more stable.
In the operational phase calibration, the phase correction coefficients are estimated using AirSWOT pixels over stationary land areas. If a high resolution digital elevation model (DEM) is available, DEM pixels over non-vegetated stationary land areas with high interferometric coherence can be used as ground control points (GCPs). The GCPs are used to estimate the phase correction needed for AirSWOT to match the GCP elevation. The estimated phase correction is added to the interferometric phase to produce a calibrated AirSWOT interferogram. If a high resolution DEM is not available, the phase correction parameters can be fitted by minimizing the elevation difference between repeat AirSWOT measurements of the same stationary surface (typically land). Essentially, pixels within these stationary surfaces are used as tie points in order to calibrate the interferometric phase of multiple AirSWOT flight lines within the overlapping area.
Either method (GCP or tie point) can be represented as a system of linear equations, solved in the AirSWOT processing chain using a weighted least squares method. However, these methods work best with the presence of bare ground areas within the study area, without significant vegetation or topography.
In vegetated areas, the AirSWOT elevation measurements are biased above ground level. The magnitude of this bias varies as a function of vegetation structure and incidence angle. This is because the microwave penetration depth within the vegetation canopy depends on incidence angle (e.g., [49][50][51]). This implies vegetated areas cannot be reliably used as either GCPs or tie points. Locations with severe topography also cannot be used, because elevation measurements from varying viewing geometries for the same location may differ as a result of foreshortening and layover [14]. The effects of foreshortening and layover in AirSWOT and SWOT data are significant due to the steep incidence angles at which these instruments operate [20,33].
We developed a new phase calibration procedure better suited to coastal wetland environments where stationary and non-vegetated land areas may be scarce or nonexistent. This new approach uses AirSWOT water surface measurements (i.e., only open water pixels) rather than land to generate a system of linear equations. The calibration approach is based on two objectives: (i) remove cross-track slopes by minimizing slopes in areas where cross-track WSS are expected to be small; and (ii) Remove absolute WSE errors by minimizing the difference between the AirSWOT WSE and in situ water level data which serve as ground control points (GCPs). The use of external elevation data for calibration of elevations estimated from radar interferometry is a standard practice [52,53], helping to account for errors in the platform position or baseline calculations, which can impact the accuracy of DEMs generated from radar interferometry. The main difference with our method is that the GCPs represent WSE rather than land elevations. Since WSE varies over time, when using the in situ WSE to calibrate AirSWOT, we interpolate the in situ WSE to the AirSWOT acquisition time using linear interpolation between the two closest measurements.
The effect of the AirSWOT phase calibration procedure on the data is illustrated in Figure 3. The original WSE measurement including the original phase drift estimate is compared to the new phase-calibrated measurements. The original WSE used a phase calibration method where the AirSWOT elevations over land were fitted to the SRTM DEM, creating a significant cross-track elevation slope, presumably due to vegetation-related biases, errors in the DEM, and a lack of land targets. To demonstrate the difference between the original calibration and our new method, we selected an area of water southwest of the Wax Lake Delta where cross-track slope should be minimal over the AirSWOT swath. In the original data the cross-track elevation slope is equal to almost 100 cm/km across the AirSWOT radar swath, an unrealistically large WSS value in this area (see Section 4.3). A higher resolution and more accurate DEM may improve those results, but the vegetation related biases would remain an issue. Here, we propose an algorithm that does not rely on a DEM. The cross-track slope is mitigated with the new calibration procedure as shown in Figure 3b.   Figure 3a shows maps of AirSWOT elevations using the original phase drift parameters (left map) and after applying our calibration method (right map). Figure 3b shows AirSWOT elevation vs. incidence angle for the water polygons (red boxes) in Figure 3a. The solid lines show the mean AirSWOT elevation at each incidence angle, while the shaded region shows ±1σ around the mean. After our calibration procedure, the anomalous cross-track slope seen in the original results has been corrected.
To preserve the true water surface slopes in the channels of the Wax Lake and Atchafalaya river systems, we fit the AirSWOT elevation profile across non-river open water surfaces to the local NAVD88 geoid (with a fitted elevation bias that varies between flight lines, to account for tidal changes). The water mask was used to select large water bodies that are not channels, as the latter may exhibit significant slopes unrelated to the geoid. A map of this calibration mask is shown in Figure 4, where black shading denotes areas that were excluded from the calibration procedure. This includes the Wax Lake and Atchafalaya rivers and deltas, as well as the areas that were classified as land in the a priori water mask. While the assumption that the water surface is locally flat with respect to the geoid in these areas may cause phase calibration errors (e.g., if the WSS is significant), we demonstrate the assumption is reasonable by showing strong agreement between calibrated AirSWOT WSEs and in situ data as shown in Section 4.
When fitting the phase correction coefficients, we fit either a zero order (constant phase correction) or first order (linear phase correction as a function of along-track distance) correction for each flight line. To fit these phase correction coefficients, the equation linking the observed AirSWOT elevation bias to the phase correction can be written as follows for a single AirSWOT image pixel: dh dφ is the conversion factor between elevation and interferometric phase for that pixel, which can be calculated using knowledge of the AirSWOT instrument characteristics and viewing and target geometry [14]. φ 0 is the zero order phase correction which we wish to solve for, and is subtracted from the measured interferometric phase to obtain the calibrated phase. ∆h represents a bias between the AirSWOT elevation and the geoid, which is held constant for each flight line. ∆h is treated as an unknown and is estimated during the calibration, as it can vary due to tide or AirSWOT GPS bias. The parameter h AirSWOT is the AirSWOT elevation estimate for the pixel, and h N AVD88 is the elevation of the NAVD88 geoid at that position. Of course, for a single image pixel this equation is underdetermined, as it has two unknowns-φ 0 and ∆h. However, since φ 0 and ∆h are constant for each AirSWOT flight line, the system of linear Equation (1) (as repeated for every AirSWOT image pixel not in the masked area of Figure 4) can be solved using least squares. Similarly to the zero order case, the case of a first order phase correction can be formulated as follows: φ 1 is the along-track slope of the differential phase correction, and s is the along-track coordinate of this image pixel in the SCH (S: along-track distance, C: cross-track distance, and H: height) coordinate system [54]. Again, while we show here the equation for a single image pixel to illustrate how the calibration is formulated, a large system of linear equations (the equation above is repeated for many AirSWOT image pixels) is generated to solve for the calibration phase. While removing the cross-track elevation slopes over the non-river areas is one objective of the calibration procedure, the other objective is to correct the AirSWOT elevation using in situ water level stations as GCPs. This approach further refines the absolute WSE calibration. While removing the cross-track elevation slope provides relative WSE calibration, incorporating in situ data provides absolute WSE calibration.
In the case of correcting the AirSWOT data with the in situ GCPs for a single AirSWOT pixel, the equations for a zero order and first order phase correction take the form: h GCP is the WSE of the GCP (interpolated at the AirSWOT acquisition time), while ∆h GCP is a fitted bias parameter, which represents the estimated WSE bias between the AirSWOT data and the absolute reference water level. River channels are included and the calibration mask is not used for absolute calibration using the GCPs to calibrate the data. However, we do mask out land areas using the a priori water mask. Around each GCP, we include pixels that cover a square area of 0.5 km 2 (a window size of 141 by 141 pixels, with pixel spacing of 5 m). If we combine the systems of linear equations for the cross-track ocean slope correction and GCP corrections, we can create a larger system of linear equations that can be used to simultaneously estimate the phase correction parameters for multiple flight lines. In this paper, we calibrated all flight lines for each day of the AirSWOT campaign as one system of linear equations. One calibration was performed for the 8 May data, one for the 9 May data, and one for the 11 May data.
Combining the above equations, the system of linear equations that is solved in order to estimate the phase drift calibration parameters, in the simplest sense, can be written as follows: If we expand A, x (the unknowns), and b (the elevation residuals) to their actual forms, we can formulate the phase drift correction estimation problem as shown in Equation (6): The phase correction and elevation bias parameters in x are then solved for using weighted least squares. The subscripted letters represent different flight lines (e.g., φ 0a represents the zero order phase correction coefficient for line a). In the example shown in Equation (6), we illustrate how the phase calibration would be formulated for a case where there are three AirSWOT flight lines (a, b, and c) and two GCPs (GCP1 and GCP2). The elevation bias parameters with respect to the NAVD88 geoid (∆h a , ∆h b , and ∆h c ) vary between flight lines (because each flight line was collected at a different time, and covers a different area). However, the absolute elevation bias with respect to the GCPs is the same for all flight lines, as we wish to use the GCPs as an absolute water level reference which we use to fit the AirSWOT data.
Each row of the matrix represents a single AirSWOT image pixel. The vector w weights each row of the matrix equation, and is equal to 1 h is the estimated AirSWOT elevation variance due to phase noise. This is calculated by the AirSWOT processor using the Cramer-Rao lower bound of the interferometric phase variance, a widely known estimate of interferometric phase accuracy in radar interferometry [14]. The estimated phase variance is scaled to a elevation variance using the value of dh dφ for each pixel.
We exclude pixels from the calibration if their elevation uncertainty σ h is greater than 3 m. These pixels tend to exhibit low interferometric coherence and are thus prone to have large uncertainty.
Our calibration includes two independent objectives: (i) the cross-track anomalous slope removal, and (ii) the GCP absolute water level calibration. However, there are a much greater number of non-river pixels than there are pixels near the water level stations used for the absolute GCP calibration, since we use only two stations (the USGS Calumet and Caillou Bay stations). Because of this difference in the number of pixels (rows in the matrix), the GCP calibration would have little effect on the results if it were weighted equally to the cross-track slope removal objective. To balance the two objectives, we multiply the weight of each of the pixels used in the GCP calibration by 100, as this made the total weight of the pixels used for the GCP calibration roughly equal to the total weight of the pixels used in the cross-track slope removal. After applying the phase calibration, we subtract the elevation bias between the AirSWOT elevations and the USGS Calumet in situ data. This bias is estimated and removed using the AirSWOT flight line over the Wax Lake Outlet on each day of the campaign. This ensures that the AirSWOT WSE closely matches the in situ water level at the Calumet station, using Calumet as the absolute elevation calibration GCP for our AirSWOT data.
The phase biases estimated with our new method for each flight line on each day of the AirSWOT campaign are shown in Figure 5. The phase biases are observed to quickly change at the beginning of each flight, becoming more stable over time. Because of this increasing stability over time, we used a first order phase correction for the flight lines in the first hour of each flight, and zero order phase corrections for the flight lines after the first hour of flight. As can be seen, the 8 May and 9 May flights contained 15 flight lines each. The 11 May flight contained only 9 flight lines, and covered a smaller coverage area over the wetlands than the 8 May and 9 May flights. The phase calibration coefficients are assumed to vary at least in part based on temperature differentials within the radar hardware. This explains why the phase calibration stabilizes over time, because we also expect the temperature to stabilize as the flight time increases. One flight line collected on 9 May is characterized by an outlier phase drift. This flight line was located at the northern edge of the study area, away from any of the in situ water level data or targets of interest. It also did not contain a significant number of non-river pixels to use in the cross-track slope removal part of the calibration. We therefore conclude that the phase calibration for this flight line was ill-posed, and we therefore excluded it from further analysis.
To demonstrate the accuracy of the new phase calibration approach, we validate the elevation accuracy of the AirSWOT data using independent in situ water level data as discussed in the following sections.

Filtering and Averaging of AirSWOT Elevation Data
The AirSWOT elevations at the pixel level require filtering and spatial averaging to improve precision and be comparable to in situ water level measurements. We estimate WSE by averaging within a square window of approximately 0.5 km 2 around each water level station (141 by 141 pixels with 5 m pixel spacing). Before averaging, we perform a number of filtering steps in order to remove pixels that would bias the average. This is particularly important because the AirSWOT elevation distribution tends to be right-tailed, such that elevation outliers will tend to positively bias the WSE estimate.
The first filtering step is to filter out land pixels, and pixels within a 10 m buffer of land. This buffer prevents inclusion of residual land pixels from water mask misregistration, changes in water level between the water mask and AirSWOT acquisitions, or other potential errors in the water mask.
Next, we identify and filter elevation outliers in two steps. First, AirSWOT pixels with elevations >3 m or <−3 m compared to the NAVD88 geoid are rejected as these values are well outside the bounds of the expected water levels within this area.
Second, we remove additional outliers using a median absolute deviation (MAD) filter [55]. The MAD filter is known to outperform the use of the mean and standard deviation for outlier removal, particularly for skewed distributions [55], such as the AirSWOT elevations. The MAD filter effectively removes additional outliers that are within the bounds of the ±3 m threshold. Its implementation computes the median elevation of the remaining AirSWOT water pixels h median , and splits the data into two sub-distributions to the left and right of the median-i.e., less than or greater than the median value. Pixels equal to the median value are included in both distributions. We then take the median value of the left and right distributions, which we call h le f tmedian and h rightmedian . A pixel is determined as an outlier using a score computed on these parameters for the left and right pixels: and h score = 0.6745 |h − h median | h rightmedian (8) That is, the absolute deviation of each pixel around the population median is divided by the median of that side's distribution, then multiplied by 0.6745. The value of 0.6745 scales the MAD to be equivalent to the standard deviation of normally distributed data [55], such that for normally distributed data, h score is equivalent to the number of standard deviations from the mean. We filter all pixels with h score > 2, removing them from the WSE estimation. Outlier filtering has also been performed by other AirSWOT studies [27,28], though using slightly different methodology.
The local average and standard deviation σ (as a metric for uncertainty) around each station is computed over the remaining pixels within the averaging window. This average AirSWOT WSE is compared to the in situ water level data. In the results we only show AirSWOT WSE for stations where the averaging window contained at least 1500 valid AirSWOT pixels after filtering (equivalent to 7.5% of the area within the window, or about 0.0375 km 2 ). We chose the value of 1500 by experimenting with different minimum pixel counts, but found that lower values resulted in larger WSE errors compared to the in situ data. For example, using a minimum pixel count of 500 instead of 1500 increased the WSE RMSE by 7 cm.
An example showing how the AirSWOT data is filtered, and the importance of these various filtering steps, is shown in Figure 6. Figure 6a shows the elevation of the AirSWOT pixels within the averaging window before (left) and after (right) filtering, while Figure 6b shows the distribution of the pixel-level AirSWOT elevations before filtering (in orange), after land and land buffer filtering (in green), and after outlier filtering (in blue). Bridges and power lines crossing the Wax Lake Outlet river channel are clearly seen in the data, and hinder accurate estimation of WSE around the Calumet station. These were successfully removed by our filtering method. In the right image of Figure 6a, the brown and purple pixels were removed by the land and 10 m land buffer masks, respectively. The orange and red pixels were removed by the fixed ± 3 m threshold and MAD outlier filtering, respectively. The fact that most of the red and orange pixels are near the edges of the land mask suggests they are affected by layover from the nearby land areas which is consistent with the AirSWOT look direction (toward the top left of the image). The layover results from several targets being at the same distance from the radar instrument, effectively combining vegetation and water surface measurements into the same radar pixel [14].   Figure 6a shows the averaging window before (left) and after (right) filtering. The left image shows the AirSWOT elevation of each pixel within the square averaging window for the Calumet station using an AirSWOT flight line on 9 May 2015. The right image shows the pixels removed by the filter in various colors as listed in the legend. The black circle at the center of each image shows the location of the USGS Calumet station. Figure 6b shows the AirSWOT elevation distribution before filtering (top), after the land and land buffer filtering (middle), and after all filtering has been performed (bottom). The dashed lines show the mean elevation of each distribution, which get closer to the in situ water level (shown as the dashed black line) as each filtering step is performed. Figure 6b, the sample mean of the selected AirSWOT pixels converges to within 7 cm of the true water level. The mean elevation of each distribution is shown as a dashed vertical line, with the in situ water level shown as the dashed black line. The mean elevation after applying the land and land buffer masks, without performing the additional filtering of outliers, has a bias of 50 cm. Due to the noise level of the AirSWOT elevation data, and the right-skewed nature of the unfiltered elevation distributions, careful filtering of the data is essential to avoid large positive elevation biases in the final results.

As shown in
We also calculate an estimate of WSE uncertainty within the averaging window. We include two sources of uncertainty: (i) random error which we compute as the standard error of the pixel-level AirSWOT elevations within the averaging window (standard deviation divided by the square root of the number of samples); and (ii) uncertainty in the vertical datum conversion from WGS84 to NAVD88 in this area, which can be calculated by NOAA VDatum [43], and is equal to 7.3 cm on average within the study area. For the sake of simplifying the calculation, we assume that these uncertainties are statistically independent and normally distributed. As discussed above, the AirSWOT elevation distribution is skewed before filtering, however, after filtering the AirSWOT elevations should be non-skewed, such that this simplification can be justified as a loose approximation. Under this assumption, the estimated AirSWOT WSE standard deviation is equal to: σ is the estimated standard deviation of the AirSWOT WSE, σ err is the standard error of the AirSWOT pixels within the averaging window, and σ datum is the standard deviation of the vertical datum conversion from WGS84 to NAVD88. This estimated AirSWOT elevation uncertainty is shown in Section 4 using error bars to display ±1σ around the AirSWOT WSE estimates.

Water Surface Elevation
Validation of AirSWOT WSE estimates with all in situ water level measurements (Section 3.2) is presented in Figure 7 as a scatter plot including all three days of the AirSWOT campaign. The source of in situ data is identified with different colors. We report the overall empirical uncertainty of the AirSWOT water levels in the form of mean absolute error (MAE) and root mean square error (RMSE) statistics. The MAE was calculated as follows: where E is the expectation operator (sample mean in this case), | | is the absolute value (applied to each element of the vector individually), h AirSWOT is a vector of the estimated AirSWOT WSEs at each water level station, and h Station is a vector of in situ water levels measured at each station (interpolated to the AirSWOT acquisition time using linear interpolation between the two closest measurements). The RMSE was calculated as follows: where the exponent is applied to each element of the vector individually.
The MAE for all validation stations was 13 cm (RMSE of 17 cm), with the largest errors observed for some of the CRMS stations. The uncertainty calculated using only non-CRMS validation stations is 4 and 5 cm for MAE and RMSE, respectively. We determined the largest errors occur at higher AirSWOT incidence angles (Figure 8). Above incidence angles of about 15 • , errors increase for the CRMS and calibration stations. The non-CRMS stations show smaller errors than the CRMS stations, and also appear to be less affected by incidence angle, though one calibration station at very far range does show a bias of 18 cm. One possible explanation for the correlation between elevation error and incidence angle is that the CRMS stations tend to be in areas with more vegetation, particularly emergent vegetation. Vegetation-related InSAR elevation biases are greater in magnitude at higher incidence angles [50,51]. CRMS stations also tend to be in areas with thinner river channels, which can make it more difficult to filter out areas contaminated by layover. Since there are fewer water pixels in the averaging window, any localized errors (such as layover) will tend to have a greater impact on the mean elevation of the pixels within the window, causing the MAD outlier filter to be less effective. In addition, errors in the AirSWOT phase calibration will result in larger elevation errors at far range than at near range, as dh dφ increases with incidence angle [14]. We found that limiting our analysis of CRMS stations to those located at AirSWOT incidence angles of less than 15 • significantly decreased the uncertainty as shown in Figure 9. The MAE decreased to 10 cm, with RMSE of 12 cm. This is only 1-2 cm larger than the uncertainty of the calibration stations. A linear fit applied to all station types has equation y = 0.99x + 0.06 with r 2 = 0.89, showing strong agreement between AirSWOT and in situ WSE. We did not restrict the incidence angle of the non-CRMS validation stations because these stations showed less of an elevation error vs. incidence angle trend, and the non-CRMS stations are more limited in number than the CRMS stations. The WSE uncertainty observed in this study is within the range of those reported in previous studies using AirSWOT. In two previous studies of the Tanana River in Alaska, Altenau et al. [27,28] reported RMSE values of 9 cm and 11.8 cm. However, we note that we use a slightly smaller window size for our spatial averaging than in those studies. We use a 0.5 km 2 window, while Altenau et al. use 1 km 2 areas. Tuozzolo et al. [30] reported RMSE of 11.6 cm compared to in situ data over the Willamette River, Oregon, when using a 1 km shifting window to perform the averaging. Pitcher et al. [29] tested window sizes of both 1 km 2 and 0.0625 km 2 for a study in the Yukon Flats Basin in Alaska, and reported RMSE values of 8 and 15 cm, respectively, for rivers only. For lakes, Pitcher et al. reported RMSE of 21 cm.
Overall, our results show a high level of agreement between the AirSWOT WSE estimates and the in situ data. The similarity in the error levels of the calibration and validation stations suggest our calibration procedure is robust and not overfitted, as it produces accurate AirSWOT WSE estimates for both the calibration and validation water level data.

Measuring Water Level Change
We assess the consistency of AirSWOT WSE measurements over time by looking at the relative change in WSE between repeat AirSWOT acquisitions. The results are shown in Figure 10. Change in WSE was calculated for all in situ water level stations observed at least twice. As in Figure 9, we excluded CRMS stations observed at AirSWOT incidence angles ≥15 • . The MAE of the WSE change estimates is slightly greater than the instantaneous WSE measurements-MAE of 11 cm and RMSE of 14 cm. The linear fit including all stations was found to be y = 1.05x − 0.04 with r 2 = 0.95, showing good agreement between AirSWOT and in situ.

Water Surface Slope in the Wax Lake Outlet
AirSWOT can generate spatially continuous maps of WSE and WSS in a very short time period. Figure 11 shows a WSE river profile along the Wax Lake Outlet for an AirSWOT flight line acquired on 9 May 2015. The profile shows AirSWOT WSE as a function of distance along the channel, with the in situ water level data shown for comparison. We calculate the WSE profile by first creating a user-defined profile line along the Wax Lake Outlet, which is shown as the red line in Figure 11b. For each AirSWOT pixel in the area we reproject the data from the SCH coordinate system [54] to an along-channel and cross-channel coordinate system using the AirSWOT geolocation information and the manually drawn profile line. To remove AirSWOT pixels far away from the river channel, we discard pixels with cross-channel coordinates less than -700 m or greater than 900 m (distances from the profile line). While these bounds are wider than the Wax Lake Outlet itself, they were chosen to include a large portion of the river delta. The bounds are asymmetrical around zero because the manually drawn profile line isn't located at the center of the river-it was drawn along the western riverbank of the Wax Lake Outlet.
To calculate the AirSWOT WSE along the profile, we performed a moving average, again removing outliers (as in Section 3.2), over a window length of 1 km in the along-channel direction, generating a WSE sample every 50 m along the profile. Again, water pixels were filtered using the land mask and 10 m buffer, followed by fixed thresholding and MAD filtering [55]. When taking the moving average, we also calculate the standard error of the AirSWOT pixels within the moving window and estimate the AirSWOT WSE uncertainty as described at the end of Section 3.2. The estimated AirSWOT WSE uncertainty is shown as the shaded blue region in the profile (±1σ around the AirSWOT WSE estimate). After calculating the moving average, we perform further smoothing of the AirSWOT profile using a first order Savitzky-Golay filter with a window size of 2 km [56]. We use this additional smoothing to reduce the effects of noise and also to calculate profiles of WSS vs. along-channel distance (see below). We perform the moving average and Savitzky-Golay filter separately because the Savitzky-Golay filter must be applied to regularly sampled data, and because performing these steps in succession performs more smoothing than if the Savitzky-Golay filter was applied directly to the AirSWOT pixel elevations. In the reference image, the land mask (in black) is superimposed on the AirSWOT elevation map for 9 May. The Wax Lake Outlet center line used to generate the profile is shown in red. The dashes along the center line correspond to the vertical dashed black lines in the profile, and represent major inflow and outflow locations along the channel. In the profile, the solid blue line shows the estimated AirSWOT WSE along the profile, while the blue shaded region shows the estimated AirSWOT WSE uncertainty (±1σ). The circular markers denote in situ water levels from the station data.
The dashed vertical black lines in the profile (Figure 11a) show the main confluence locations along the channel, which correspond to the dashes along the red line in the reference image below the profile (Figure 11b). The first confluence is the Intracoastal Waterway (ICW). Upstream of the ICW, the Wax Lake Outlet has a significantly steeper WSS than downstream of the ICW.
The profile shows the in situ water levels for the USGS Calumet, CRMS 0479, and LSU Mike Island stations. At the downstream end of the profile, the AirSWOT WSE is slightly overestimated compared to the station data, though still within the estimated AirSWOT WSE uncertainty.
On 8 May and 11 May, AirSWOT WSE profiles along the Wax Lake Outlet have only partial coverage due to low radar backscatter. This phenomenon has been observed in several campaigns and is generally referred to as "dark water". It results from specular reflection of microwaves over a smooth water surface often encountered in very low wind conditions. The σ 0 of "dark water" is close to or below the AirSWOT instrument noise level, resulting in very low interferometric coherence, and thus noisy interferometric phase. Nonetheless, partial profiles were extracted for 8 May and 11 May ( Figure 12) and suffice to estimate large scale WSS.  (Figure 12b) showing the effects of "dark water" on AirSWOT WSE. "Dark water" results in void data and potential anomalies in the small, disconnected regions between the void data.
We note that in some of the areas where the profile could be estimated, erroneous WSE measurements (e.g., slopes in the wrong direction) are found at the edges of disconnected segments of the WSE profiles. A possible explanation is that while the water was specular, layover signals caused by vegetation along the river bank were interpreted as water surface signal. Another explanation is that when the amount of valid water data is small, the profile is more sensitive to errors in the water mask which allow land pixels to be incorrectly included in the profile. Despite the discontinuities in the WSE profile for two days of the campaign, we can use the AirSWOT measurements to calculate the overall WSS along the river channel between the USGS Calumet and LSU Mike Island stations.
A comparison of in situ vs. AirSWOT-derived WSS estimates for the Wax Lake Outlet are shown in Table 3. To compute these slopes, we simply took the rise over run for the AirSWOT and in situ data at the Calumet and Mike Island stations. The largest error was for the 9 May flight, with absolute WSS error of 0.29 cm/km, corresponding to a percent error of 7.7%. The other two flights had smaller slope errors: −0.041 cm/km (−1.0%) on 8 May, and 0.17 cm/km (3.8%) on 11 May. These WSS errors are significantly smaller than those reported in previous studies using AirSWOT to measure WSS [27][28][29][30], which generally found WSS errors on the order of 1.0 cm/km. However, we note that while those previous studies calculated WSS over 10 km reaches, here we calculate the slope between the USGS station at Calumet and the LSU Mike Island water level gauge, which are separated by almost 24 km. Unfortunately, we cannot accurately validate WSS over 10 km reaches of the Wax Lake Outlet due to the sparse coverage of the available in situ water level data in this location, as well as the void data caused by specular water on the 8 May and 11 May flights.
In general, we expect longer river reach distances to yield smaller WSS errors, as has been shown by Altenau et al. [28]. This is the case because as the river reach distance increases, the denominator of the slope equation (the reach distance) increases, but the error in the numerator of the slope equation (the change in water level along the reach) remains the same. Another explanation for these small WSS errors is that the Wax Lake Outlet is relatively straight, such that the AirSWOT incidence angle range was small between Calumet and Mike Island. This reduces the effect of any incidence angle-dependent errors, such as elevation biases due to residual errors in the phase calibration. The fact that our calibration procedure uses in situ data may also contribute to our better results.
One main advantage of AirSWOT or other remote sensing techniques compared to in situ data is the spatial continuity of the WSE and WSS estimates. In addition to deriving an average WSS between stations, AirSWOT can produce a profile of WSE and WSS along a river channel, in particular when there are no "dark water" issues. A WSS profile for the 9 May data is shown in Figure 13. The blue line shows WSS estimated using a first order Savitzky-Golay filter with a 2 km moving window, the same window size used to produce the WSE profiles. We also produced a WSS profile using a 10 km moving window, to match the 10 km reach length used in other AirSWOT studies [27][28][29][30]. The WSS profiles do not extend to the full length of the WSE profiles due to the length of the Savitzky-Golay filter.
AirSWOT-derived WSS of the Wax Lake Outlet was steepest near the USGS Calumet station, reaching a value of −14 cm/km before the slope decreased in magnitude as water approached the ICW (denoted by the first dashed black line). Downstream of the ICW, WSS mostly varied between −4 to −1 cm/km. It is not possible to validate these fine resolution WSS estimates without a dense set of in situ water level measurements. A future study should plan for a spatially dense set of in situ measurements over the primary river channels, similar to other studies [27][28][29][30].

Conclusions
We have investigated and demonstrated the use of AirSWOT to map WSE and WSS in coastal marsh environments. Our results show AirSWOT's capabilities as a science instrument to support the upcoming Delta-X campaign and SWOT calibration/validation efforts. SWOT will provide a global perspective of WSE and WSS when it launches in 2022, while the airborne platform of AirSWOT enables focused studies at fine spatial and temporal resolution. This makes AirSWOT a particularly valuable instrument in highly dynamic and spatially complex coastal environments.
We implemented a data-driven method to calibrate the AirSWOT interferometric phase, which estimates and then removes time-varying differential phase delays between channels in the radar hardware. Our method uses non-vegetated water surfaces to remove instrument-related cross-track elevation slope and in situ water level measurements to remove remaining biases. Our method can be applied to future AirSWOT campaigns in coastal deltas and estuaries (e.g., the Delta-X campaign). While other AirSWOT studies [27][28][29][30] have focused on science results, this study introduces the details of AirSWOT phase calibration.
After applying the phase calibration method, we validated AirSWOT WSE in open water and found RMS errors of 12 cm over areas less than 0.5 km 2 . The SWOT science requirement for WSE accuracy is 10 cm RMSE averaged over a water area of 1 km 2 . While our RMSE is slightly larger, we averaged over a much smaller area to produce this error level. We only included results in our analysis if they were averaged over at least 1500 AirSWOT pixels, which covers a water area of only 0.0375 km 2 . Our error levels are similar to that reported in previous AirSWOT studies [27][28][29][30].
We used AirSWOT to estimate profiles of WSE and WSS along the Wax Lake Outlet river channel, with slope error of less than 0.3 cm/km on all three days of the campaign, less than 8% of in situ WSS. AirSWOT can estimate WSS with very high accuracy. The slope errors we reported are less than in previous studies [27][28][29][30], but we noted that the Wax Lake Outlet is a very straight river channel, reducing the potential effect of incidence angle-dependent errors (i.e., phase calibration errors). This is the first study using AirSWOT to estimate WSE and WSS in coastal deltas and wetlands. Despite the challenges involved when processing AirSWOT data in these complex and spatially heterogeneous environments, our results demonstrate that AirSWOT provides accurate estimates of open water WSE and WSS in coastal deltas. The methods and results from this study will be used to support the future Delta-X and SWOT calibration campaigns.
Author Contributions: M.S. and E.R. conceptualized the study. M.S. supervised the study and secured funding. M.D. developed the methodology for AirSWOT phase calibration in coastal wetlands and for filtering and spatial averaging of the data, and implemented these methods. M.D. performed investigation and analysis of the remote sensing and in situ datasets. X.W. and A.C. worked on the AirSWOT processing software and assisted with the processing. M.D. prepared the original draft of the paper. All authors (M.D., M.S., E.R., X.W., A.C., and T.P.) contributed to the revision and editing of subsequent drafts.

Funding:
The research described in this paper was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. This project was supported, in part, by NASA's Physical Oceanography as well as the Terrestrial Ecology programs.