Impacts of Spaceborne GNSS-R System Observations on Storm Surge Prediction: A Hurricane Harvey Case Study Using CYGNSS Full DDM Downlinks

This work assesses the impact of the Cyclone Global Navigation Satellite System (CYGNSS) (full delay-Doppler map) informed wind fields on storm surge simulations using the constellation's storm observations. In order to assess the impact of the CYGNSS-enhanced wind fields, storm surge simulations are performed with the ADvanced CIRCulation (ADCIRC) model and validation studies are performed with high water mark data provided by the U.S. Geological Survey. To provide context for the CYGNSS-based results, comparisons to storm surge predictions using Modern-Era Retrospective analysis for Research and Applications, Version 2,wind fields are also performed using the example of hurricane Harvey. In this initial assessment, it is observed that augmenting existing wind estimates with information provided by Global Navigation Satellite Systems Reflectometry systems, such as CYGNSS, has the potential of improving surge predictions relative to existing sources of wind information.

waves, tides, and rainfall.Numerical models, such assea, lake, and overland surges from hurricanes [2] and ADvanced CIRCulation (ADCIRC) [3] can be applied to predict and understand the storm surge caused by hurricanes.
Storm surge models require input information on atmospheric forcing that can be acquired from various types of datasets.Data-assimilated wind analysis products, such as the real-time hurricane wind analysis system (HWind) [4] combine various observations of wind velocity from land-, air-, and spacebased platforms [5].Alternatively, atmospheric models, such as the weather research and forecasting (WRF) Model [7] and the regional atmospheric modeling system [6] are based on the solution of physics-based governing equations.Such models are often further coupled with ocean circulation/wave models to provide a better representation of wind fields as in the unified wave interface-coupled model [8] and the Northeast Coastal Ocean Forecasting System.In contrast, reanalysis products are produced by combining the physics-based and data-assimilation approaches as in the HWind and interactive objective kinematic analysis system [9] of Oceanweather Inc. [10].
A preferred "parametric model" approach describes wind and pressure fields using analytic expressions (e.g., [11], [12], and [13]) based on a small number of storm parameters, such as maximum wind speed, location of storm center, and radius to the maximum wind speed, simplifying the required input information.The generalized asymmetric Holland model (GAHM) [13] is one such model that is widely used for storm surge simulations.The storm parameters that inform the model are typically obtained from the National Hurricane Center's advisory (for forecasting) and best-track (for hindcasting) data.Recent studies [14], [15], [16], [17] have demonstrated the efficacy of the parametric model approach in storm surge modeling.
Previous parametric-model-based storm surge studies have emphasized the use of best-track data for the required storm parameters.Best-track data are obtained primarily from airborne measurements by the U.S. Air Force and National Oceanic and Atmospheric Administration, with ship reports, surface observations from land stations, and data buoys also used [18], [19] along with space-based measurements to a lesser extent.
Because storms develop and undergo rapid intensification over the open ocean where in-situ and airborne observations are less frequent, an increased use of spaceborne measurements in storm surge predictions is well motivated.NASA's 2016 Earth Venture Cyclone Global Navigation Satellite System (CYGNSS) Mission currently operates a seven satellite constellation [20], [21] performing ocean wind speed measurements using a mode of microwave remote sensing known as Global Navigation Satellite System Reflectometry (GNSS-R).This article reports results from a study of the impact of the CYGNSS-enhanced wind fields on storm surge simulations.ADCIRC storm surge simulations are compared for input wind fields including or excluding CYGNSS information, and validation studies are performed with high water mark (HWM) data provided by the U.S. Geological Survey (USGS).

II. ADCIRC FORMULATION
Storm surge models are typically based on the depth-averaged shallow water equations (SWEs), which are derived from the Navier-Stokes equations following vertical integration and an assumption of hydrostatic pressure.The depth-averaged SWEs consist of the continuity equation and the momentum equation where ζ is the free surface water elevation from the geoid, h b is the bathymetric depth relative to the geoid, H = ζ + h b is the total depth of water column, q ≡ (Hu, Hv) where u and v are the depth-averaged velocities in the x and y directions, respectively, τ b and τ s are the bottom and free-surface stresses, respectively, ρ is the density of water, f c is the Coriolis parameter, g is the gravitational acceleration, ε is the eddy viscosity coefficient, and F represents any additional body and surface forces.
In ADCIRC, the bottom and free-surface stresses are computed by standard quadratic drag laws where C b is a bottom friction coefficient, C w is a wind drag coefficient, ρ air is the air density, and w = (w x , w y ) is the wind velocity.The bottom friction coefficient can be given by the Manning or Chezy coefficient, and the wind drag coefficient is defined by Garratt's drag formula [22] C w = (.75 + .06w ) × 10 −3 . ( ADCIRC solves a "wave continuity equation" introduced by Gray and Lynch [23] and extended to a generalized wave continuity equation (GWCE) by Kinnmark [24].The wave continuity equation and GWCE are reformulations of the continuity equation, obtained by combining (1) and (2).ADCIRC solves the GWCE with a continuous Garlerkin approach that provides monotonic dispersion relationships and eliminates spurious oscillatory solutions.As in [37], the results to be shown also couple ADCIRC with the Simulating WAves Nearshore (SWAN) model that describes short-crested waves in coastal regions [36].The ADCIRC and SWAN models are coupled in two-way communication: ADCIRC computes wind velocities, water levels, and currents and passes them to SWAN, and then SWAN computes wave radiation stress gradients that are added to ADCIRC's free surface stress term in (2).ADCIRC+SWAN performs storm surge simulations that include tides and waves to improve prediction accuracy [38].

III. MATCHED FILTER WIND RETRIEVALS
For the high wind speeds of interest in hurricane studies, the baseline wind speed retrievals reported by CYGNSS in its standard observation mode exhibit increased errors.Due to this limitation, a "matched filter" wind speed retrieval approach that is more applicable for hurricane observation scenarios [25], [26] is adopted.The matched filter retrieval reports maximum sustained surface winds (SSW) V max by comparing CYGNSS delay-Doppler map (DDM) measurements in the "full DDM" CYGNSS special observing mode over a track of multiple specular points.Measured "full DDM" values are compared with those predicted by an end-to-end simulator forward model [27], [28], [29], [30] of the same.The track of observations is determined using a series of CYGNSS/GPS Earth Centric Earth Fixed position and velocity vectors which are used to model forward signal propagation, scattering off the ocean surface, and mapping of surface scatter into delay-Doppler space for a parametric cyclone model [31] of varying V max .An example comparing simulated and measured CYGNSS DDMs over comparable surface conditions is depicted in Fig. 1 highlighting the ability to accurately describe DDM power distributions accurately using the forward modeling architecture.
To retrieve a maximum wind, both the measured M p (τ, f D ) and simulated S p (τ, f D ; V max ) DDMs are initially normalized (i.e., all DDM reported values are divided by the maximum amplitude within the DDM) and the correlations across a set of P specular points is computed over an SSW range of 15 ≤ V max ≤ 80 m/s where Fig. 2 depicts R 2 values as a function of the parametric storm V max for an example track.While the variability in storm structure with V max produces only modest changes in along track correlations, the results typically show a clear peak that identifies the "retrieved" V max value.
While this "track-based" retrieval approach is significantly more complicated than the "point-wise" standard Level-2 wind speed retrieval, it offers multiple advantages.Its reliance on "DDM overall shape," as opposed to amplitude, over a track reduces the impact of CYGNSS DDM amplitude calibration errors.The extended delay range of full DDM observations and comparison with a forward model of the same over an extended storm overpass enables information on storm structure to be incorporated into the retrieval.The requirement for full DDM mode measurements however reduces the coverage available, as CYGNSS receivers are only occasionally commanded to operate in this mode.

A. CYGNSS Maximum Wind Retrievals
Hurricane Harvey emerged from the African west coast as a tropical wave on DOY 224, 2017, and over its ≈30 day life cycle underwent a sustained SSW increase to a maximum of 62 m/s and a category four characterization on the Saffir-Simpson scale on DOY 238, 2017.The storm's landfall in Texas on DOY 237 and in Louisiana on DOY 242 induced peak storm surge levels of 3.6 m and caused approximately $125 billion of damage [19].
Numerous CYGNSS full DDM mode [32] collections targeting hurricane Harvey were performed beginning from DOY 231, 2017.Fig. 3(a) illustrates CYGNSS specular points for which full DDM measurements targeting Harvey were performed on DOY 232 overlaid on best-track storm center data for the entire storm history, with positions on DOY 232 marked with markers.full DDM measurements occur as tracks of individual specular points that can span in excess of 1000 km.Because the exact location of the storm center is not known at the time the full DDM collection is scheduled, many full DDM measurements occur several hundred km from the storm center, excluding their use in the matched filter retrieval process.A CYGNSS/storm center time/space colocation criterion was enforced to identify measurements having significant overlap with the storm's major features.Hurricane Harvey storm center locations were initially linearly interpolated within the 6 h best-track reporting interval, and locations occurring within ±15 min of a CYGNSS tracks were further used to apply a 100 km colocation criteria at minimum separation.For each track satisfying these criteria, 51 specular points (each within approximately 150 km of the closest location) centered about the measurement closest to the storm center were used in the "track-based" retrieval.The retained tracks are depicted in Fig. 3(b) over the DOY 231-238 period when full DDM downlinks were available.
The ADCIRC simulations to be performed require input storm parametric wind fields at ≈ 6 h intervals over the storm's life cycle.Because CYGNSS full DDM mode retrievals are not available on this time scale, an alternative approach based on augmenting Best track maximum wind estimates with matchedfilter V max retrievals is explored.Consider a best track maximum wind time series V BT max (t) and a CYGNSS maximum wind retrieval at time t 1 V C max (t 1 ); an "augmented" time series V A max (t) is formed by updating all maximum wind best-track time estimates following time t 1 and prior to a subsequent CYGNSS measurement at time t 2 by a factor proportional to the ratio of CYGNSS/best track maximum winds as follows: The process is repeated across all available CYGNSS measurements until a complete V max record is formed in preparation for use as part of ADCIRC's storm surge simulations.While it is acknowledged that the resulting set of wind fields are not purely representative of CYGNSS observed winds given the rescaling of best-track time series as part of the process, related improvements or degradation in storm surge results relative to in situ gauge data are nonetheless expected to be indicative of  the potential or lack thereof for improving surge level prediction in an operational capacity as future GNSS-R receivers and extended full DDM mode coverage come online, further shortening revisit time and expanding coverage such that more independent data records may be constructed.
The resulting CYGNSS-enhanced V max time-series using the eight available CYGNSS retrievals is depicted in Fig. 4; times of the associated CYGNSS measurements are indicated by the filled markers.The results show hurricane Harvey's development from a tropical depression to a Cat. 4 hurricane.An initial visual examination of the original best-track and CYGNSS augmented time series makes clear the complementary nature of both sets of maximum wind estimates while a more pronounced difference is observed on DOY 237 where the CYGNSS derived results suggest higher maximum winds of up to 62 m/s compared with the reference best track estimate of 50.29 m/s at a comparable timestamp.The resulting peak 11.71 m/s difference increases subsequent portions of the time series by a factor of 23.28%, while the average scale factor is 10.25%.In addition, the two time series show a mean V max difference of 3.50 m/s, a root mean square error (RMSE) of 3.78 m/s, an unbiased-RMSE of 3.50 m/s, and a mean retrieval error relative to the reference best track estimate of 11.88%.This relative difference is comparable to the maximum 10% targeted error levels in high wind speed retrievals identified as one of the CYGNSS mission requirements.The overall correlation of the two time series is 90.24%.

B. Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA-2) Maximum Wind Retrievals
MERRA-2 reported winds were also used to provide context for the impact of the CYGNSS derived maximum winds.In contrast to winds derived from individual spaceborne sensors, the MERRA-2 product is a ≈50 km gridded product that is reported on an hourly basis.A region spanning 250 km about the interpolated storm location was first identified and the best-track time series was then rescaled by the maximum reported wind field value within this subset of pixels.The resulting MERRA-2 time series is also depicted in Fig. 4, and shows a consistent underestimation of the storm's maximum winds (an average 80% downscaling of the best track V max ).For this case, the mean V max retrieval difference and RMSE were 10.48 m/s, unbiased-RMSE is 7.24 m/s, and mean retrieval error relative Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.to the reference best track estimates is 43.72%.The overall correlation of the two time series reduces to 18.14%.
MERRA-2 maximum winds were also highly consistent irrespective of storm development, exhibiting a standard deviation of 1.79 m/s.Both the underestimation of maximum surface winds and their consistency are attributable to the processes by which the MERRA-2 product assimilates a large number of observations for a given location and time period.As part of this process, ensuring temporal and spatial homogeneity are crucial such that available estimates undergo a process of recursive filtering [35] for every wind field cycle, which in turn significantly improves the quality of overall reported winds in the context of general ocean surface wind analysis but adversely impacts cases were true anomalies exist.

V. STORM SURGE SIMULATIONS
Storm surge simulations for hurricane Harvey were performed using the GAHM parametric storm model with V max values determined from the best-track, CYGNSS, and MERRA-2 timeseries shown in Fig. 4. It is noted that the GAHM model's ability to provide more realistic storm structure representations is in part enabled by its ingestion of a large number of ancillary parameters, besides V max , obtained from the best-track database.For more details see [26].Simulations were performed from 2 August 0000 UTC to 28 August 1200 UTC, 2017, including the time of landfall of Harvey (on San Jose Island, TX, around 26 August 0300 UTC).The simulations included spin-up periods of 21.5 days with meteorological forcing beginning 23 August 1200 UTC.The mesh used for the simulations is the so-called TX2008, which has high resolution near San Jose Island (see Fig. 5).
Maximum water elevations obtained from the AD-CIRC+SWAN simulations are shown in Fig. 6.It can be seen that CYGNSS-based wind fields result in a similar level of accuracy to that of best-track-based wind fields, whereas MERRA-2based wind fields result in much lower surge predictions.Validation studies were performed with HWMs provided by USGS, which are post-flood measurements of the highest elevation of floodwaters.They are measured on a number of points and thus can be used to estimate the overall accuracy of storm surge simulations.Although a total of 2364 HWMs were measured for Harvey, 1,136 of the measurements were located outside the TX2008 mesh, and an additional 1075 points were impacted by  The comparison further considers ADCIRC location classified as "partially wet," in which one node of the applicable "element" is never inundated (dry node) while another node is inundated (wet node) during the simulations.For such elements, peak water elevations are computed by assuming the water elevation equals the mesh bathymetry of the dry nodes.This computation allows estimation of water elevations on partially wet elements without requiring additional information.
Fig. 7 provides scatter plots of ADCIRC and measured HWMs for the three input wind fields.The results generally show underprediction for all three wind fields that are likely related to the approximate 1.53 m of hurricane Harvey's precipitation [19] that is not included in the simulations.Each comparison is summarized using three values: the slope of the best-fit line, R 2 , and RMSE, with the resulting values provided in Table I.The results shows that ADCIRC predictions using the CYGNSS-enhanced wind fields show the best performance for all three summary metrics, including an RMSE of 0.573 m, while the original best-track wind fields show reduced (but similar performance) and MERRA-2 winds result in significantly degraded results.

VI. CONCLUSION
The results of the hurricane Harvey case study of this article show that storm surge predictions can be improved by incorporating CYGNSS full DDM mode V max retrievals.The increased V max values obtained near storm landfall appear to be the primary factor driving these results.Although the results obtained are limited by the relative sparseness of full DDM retrievals for the period considered, the improved performance achieved is encouraging with regard to the future potential use of GNSS-R wind retrievals in storm surge prediction.Continued applications of the method are expected in the future are the dataset of available full DDM storm observations continues to grow.Future work will also include application of the proposed methodology to other storms having a variety of characteristics and for areas having varying bathymetric features.

Fig. 2 .
Fig. 2. Example matched filter correlation profile versus reference library maximum wind.

Fig. 3 .
Fig. 3. (a) Positions of all full DDM specular measurement made on DOY 232, 2017.(b) Subset of specular points meeting time/space colocation criteria over all full DDM downlinks made for hurricane Harvey.

Fig. 4 .
Fig. 4. Complete hurricane Harvey maximum wind time series input to ADCIRC.Points, at which rescaling of best-track winds takes place are indicated by the filled markers.

Fig. 7 .
Fig. 7. Scatter plot of USGS HWMs and peak water elevations of ADCIRC+SWAN from wind fields based on (a) best-track, (b) CYGNSS, and (c) MERRA2.Red solid line is best-fit line and back dashed lines indicate relative error with increment 10%.Red and blue circles indicate overprediction and underprediction, respectively.

TABLE I SUMMARY
TABLE OF REGRESSION ANALYSIS OF ADCIRC+SWAN SIMULATIONS RESULTS WITH THREE DIFFERENT WIND FIELDS riverine flooding not modeled by ADCIRC and/or were more than 200 km from the landfall location.The remaining 153 HWMs were used in the performance assessment.