## Abstract

On 26 August 2017, Hurricane Harvey struck the Gulf Coast as a category four cyclone depositing ~95 km^{3} of water, making it the wettest cyclone in U.S. history. Water left in Harvey’s wake should cause elastic loading and subsidence of Earth’s crust, and uplift as it drains into the ocean and evaporates. To track daily changes of transient water storage, we use Global Positioning System (GPS) measurements, finding a clear migration of subsidence (up to 21 mm) and horizontal motion (up to 4 mm) across the Gulf Coast, followed by gradual uplift over a 5-week period. Inversion of these data shows that a third of Harvey’s total stormwater was captured on land (25.7 ± 3.0 km^{3}), indicating that the rest drained rapidly into the ocean at a rate of 8.2 km^{3}/day, with the remaining stored water gradually lost over the following 5 weeks at ~1 km^{3}/day, primarily by evapotranspiration. These results indicate that GPS networks can remotely track the spatial extent and daily evolution of terrestrial water storage following transient, extreme precipitation events, with implications for improving operational flood forecasts and understanding the response of drainage systems to large influxes of water.

## INTRODUCTION

Hurricane Harvey migrated a total of ~900 km east along the Gulf Coast, first stalling for 2 days after making landfall in southwest Texas on 26 August 2017, after which it retreated offshore before making landfall a second time east of Houston on 30 August (*1*). Over a period of 7 days, Harvey produced record-breaking rainfall, with 1.54 m of precipitation measured southeast of Houston, causing extensive flooding and $125 billion in damage, an amount second only to Hurricane Katrina (*2*). Field surveys, satellite radar, and optical images provided measurements of water extent and depth, but these approaches could not continuously track the evolution of terrestrial water storage (TWS), which describes the sum of surface standing water, soil moisture, and groundwater. Over the past 15 years, changes in TWS have been inferred from the measured changes in Earth’s gravity field by the Gravity Recovery and Climate Experiment (GRACE) satellites (*2*). However, GRACE measurements are too coarse in both space (~300 km) and time (1 month) to accurately quantify TWS changes following transient and extreme hydrologic events such as hurricanes. Measurements of TWS at higher spatial and temporal resolution are critical to understanding how the hydrologic system responds to extreme precipitation events and for monitoring the continued and secondary flood hazard posed by stored water to downstream rivers and dams (*3*).

Water loading of Earth’s elastic crust produces instantaneous elastic deformation, which can be measured with millimeter-level precision using the Global Positioning System (GPS) (*4*). The elastic response of solid Earth to hydrological loading has previously been used to infer monthly to interannual changes in TWS such as monsoonal effects (*5*), monthly variations from seasonal rainfall in the Pacific Northwest (*6*), and water loss of the 2011–2015 California drought (*7*, *8*). Because of the rapid decay of displacement with distance from surface loads (*4*, *7*), GPS measurements can resolve fine-scale spatial variations in TWS (on the order of 10 to 100 km) and over wide, potentially difficult to reach regions. Where highly precise GPS positions are estimated daily, this approach provides a means to evaluate transient changes of water storage over large regions and with relatively high spatial resolution.

Here, we use daily GPS measurements of vertical and horizontal crustal motion to track the evolution of TWS as it accumulates and dissipates during and following Hurricane Harvey and study its effect on the region’s hydrologic system (see the “Processing of GPS data” section). However, resolving daily variations of TWS from transient hydrologic events is challenging, largely because the vertical component of GPS, which records the primary response of Earth’s crust to water loading, has the highest noise level (1σ = 5 mm, compared to 1σ = 2 mm for the horizontal). To resolve this, we use independent component analysis (ICA), a spatiotemporal filtering technique that allows us to extract the transient hydrologic signal and remove systematic biases (*9*). We apply this technique to GPS time series from all available 219 stations over a ~2-month period, from 1 August to 10 October, that are distributed every ~20 km across south Texas and Louisiana (Fig. 1). To further help determine the hydrologic loading signal, corrections for nontidal oceanic and atmospheric loading are also applied to the GPS time series (*10*).

## RESULTS

### GPS filtering—ICA

ICA is a form of blind source separation that seeks to separate underlying sources assumed to be statistically independent (for more details, see the “Spatiotemporal filtering—ICA” section) (*9*). Here, we use ICA to separate the hydrologic loading signal from temporally incoherent, spatially correlated noise, termed common-mode error (CME). Although the sources of CME are still debated, it is thought to result from a combination of uncertainties in GPS orbital, reference frame, and large-scale atmospheric modeling when estimating the GPS position (*11*–*13*) and typically contributes the largest error to GPS data (*14*, *15*). Choosing to decompose the data into four independent components (ICs) (informed from two stopping methods, see Methods), we identify the first IC as CME as it (i) contributes the largest variance, (ii) has a common spatial response across all stations, and (iii) exhibits temporally incoherent motions with a spectrum consistent with flicker noise (slope of −1; see Methods, text S1, and figs. S1 and S2) (*12*). For the vertical data, the hydrologic loading is found on the second and fourth ICs (fig. S3), associated with the first and second landfall of Harvey, respectively. Although the ICs are separated on a statistical basis and not a purely physical one, with the potential for mixing between them, the separation of both landfalls onto different components is not surprising, given that they occurred in spatially distinct areas and at different times. The third IC is a linear trend affecting stations around Houston that is likely related to ongoing groundwater extraction (see text S1). Selecting only the second and fourth ICs to define the spatiotemporal variations of the vertical GPS time series reduces the weighted root mean square (WRMS) by 55% (a median repeatability of the raw time series of 5.58 mm, and 2.49 mm for the filtered). Similarly, selecting two ICs each for the east and north components (see figs. S4 and S5) reduces the WRMS by 75 and 76%, respectively (a median repeatability of the raw time series of 1.68 and 1.90 mm for the east and north motions, respectively, and 0.42 and 0.45 mm for the filtered). Selecting components that capture the hydrological loading also serves to remove motions local to each GPS station, which can arise from monument instability, anthropogenic processes, and multipath scattering (*13*, *14*).

Humidity changes in the troposphere that can delay the GPS satellite-receiver signal and poroelastic uplift induced by aquifer recharge can both act to potentially bias the GPS motions from hydrologic loading. However, we note that these effects are unlikely to significantly alter the GPS motions, as the former is modeled in processing using weather model data (*16*, *17*) and estimated every 5 min, and the latter is inconsistent with surface motions predicted from poroelastic recharge as estimated from 30 wellhead measurements from aquifers across Houston (fig. S6). In addition, we find strong agreement between the filtered GPS motions and predictions from the North America Land Data Assimilation System (NLDAS) hydrologic model [a median Pearson correlation coefficient (ρ) of 0.72], the latter simulating hourly changes in water storage driven by observed precipitation (see fig. S7) (*18*). However, there is notably weaker agreement in the horizontal direction between the observed and NLDAS-predicted motion (median ρ of 0.52 and 0.27 in the north and east direction, respectively; figs. S8 and S9). This lower correlation could be the result of water loads being placed on the wrong side of a GPS station in the NLDAS model, due either to assumptions of how water is routed across the terrain via surface flow or rivers, or from NLDAS averaging water over 1/8-degree cells. Either of these would invert the predicted motions in time series, producing an anticorrelation between that observed and simulated, which can be seen in some cases as high negative ρ values in figs. S8 and S9. Although there are differences in the functional form of the horizontal positions, the first-order agreement in the overall relative amplitude and timing of position change is encouraging, as the horizontal component is largely insensitive to tropospheric delays, giving confidence that the observed deviations likely reflect changes induced by hydrological loading.

The filtered GPS data reveal a clear migration of subsidence across the Gulf Coast, followed by gradual uplift after Harvey dissipated (Fig. 2). One day after Harvey made landfall, 67 GPS stations around Houston subsided by up to 21 mm, with deviations of up to 4 mm in both horizontal directions (figs. S8 and S9). Five days later, 13 stations ~200 km to the east along the Texas-Louisiana border showed smaller amounts of subsidence (up to 7 mm). After 1 September as Harvey dissipated, GPS stations across the entire Gulf Coast showed rapidly decaying uplift and horizontal motions opposite to that which occurred during landfall, with the fastest uplift rates around Houston, and most stations returning to their pre-Harvey elevations and horizontal positions after ~5 weeks (Figs. 1 and 2 and figs. S7 to S9).

### Inverting GPS motions for TWS

From the filtered vertical GPS time series, we inverted all three components of motion for daily changes in TWS on a 0.25° grid spanning Texas and west Louisiana. At each grid node, we estimated water volume by minimizing the difference between the observed vertical and horizontal GPS motions with displacements predicted from a model of water loading that deforms a spherical, vertically layered elastic Earth structure (Fig. 2) (*4*, *19*, *20*). We regularize the solution by applying spatial and temporal smoothing, choosing values that minimize the reduction in the percent of variance the model can explain (fig. S10 and text S2). As informed by synthetic tests, these smoothing values are sufficient to recover loading sources of >50 km in wavelength and rates of change of water volume similar to that expected from Harvey (see figs. S11 and S12 and text S3). The preferred inversion result yields spatially random residuals, with the predicted result during Harvey’s occurrence reducing the WRMS of the observed motions by 54%.

Our inversion of the GPS data reveals that a third of the total water deposited by Hurricane Harvey (95 km^{3}) was captured and temporarily stored on land (Fig. 3), reaching a maximum storage capacity of 25.7 ± 3.0 km^{3} by 1 September (Fig. 4). The TWS estimates show two main regions that correspond to the two regions of observed GPS subsidence (Figs. 2 and 3). The first region developed when Harvey made initial landfall, loading an area that covers much of Houston, that grew steadily in volume to 22.33 ± 6.0 km^{3} over a 2-day period as Harvey stalled to the southwest. The second region developed 4 days later when Harvey made landfall a second time east of Houston, where TWS increased to 6.5 ± 1.1 km^{3} along the Texas-Louisiana border in a north-south pattern. We note that an increase of TWS of 5 to 6 km^{3} before Harvey from 7 to 9 August corresponds to a smaller known storm that occurred around the Houston area.

Precipitation data from National Oceanic and Atmospheric Administration (NOAA) stage IV (*20*) show a consistent spatial and temporal pattern with the accumulation of TWS—the former a flux and the latter an integrative quantity (Fig. 3)—giving confidence that the inversion result is reflecting hydrologic changes. In addition, comparison of our area-integrated TWS shows strong agreement in functional form and amplitude with the total TWS simulated from the NLDAS hydrologic model (green and black lines in Fig. 4A; ρ = 0.85; variance reduction, 89.37%) (*18*). Comparison of our TWS estimates at the site of the Barker and Addicks reservoirs, in Houston, also shows strong agreement with the measured reservoir water volume (ρ = 0.69; variance reduction, 65%; Fig. 4C). Last, groundwater wellhead measurements from 30 aquifers across Houston help confirm the GPS subsidence and uplift occurred synchronously with the addition and loss of water, respectively (Fig. 4B).

## DISCUSSION

To understand how a region’s hydrologic system responds to large influxes of water, and how it may behave in future events, we must assess the behavior of the components of the system, namely, water input (precipitation), storage (TWS), and output [surface runoff, groundwater flow, and evapotranspiration (*ET*)]. After 1 September when stored water reached a maximum capacity of 25.71 ± 3.0 km^{3} (1σ), it dissipated gradually at a rate of ~1 km^{3}/day over the following 5-week period (Figs. 3 and 4A). *ET* estimated from solving for surface heat flux with constraints from satellite and ground measurements (*21*) shows similar rates of water loss in the region, ranging from 0.5 to 1.5 km^{3}/day and averaging 1.1 km^{3}/day, indicating that it contributed to the majority of gradual loss of stored water over time. In addition to constraining the gradual loss of stormwater, we can also estimate the rapid phase of water loss in the days following Harvey’s arrival, namely, from the sum of the surface runoff and groundwater flow (*SG*) components. Following mass conservation, we can estimate daily *SG* by simply closing the water budget (*SG* = *P* − *ET* − Δ*TWS*) with daily estimates of *ET* (*21*), precipitation (*P*) (*20*), and changes in water storage (Δ*TWS*), which essentially provides a stormflow hydrograph of *SG*. Stormflow hydrographs are typically derived from streamflow discharge measurements that constrain the loss of stormflow waters via rivers but can be problematic to interpret as they include the effects of “old” preexisting water (*22*) and rely upon empirical relations between measured river level and discharge constrained during ordinary flow, which can be unreliable at extreme values during a large storm event (*23*). The GPS-derived *SG* component, which is not subject to these issues (cyan line in Fig. 4A), indicates that water was lost initially at a rate of ~8.2 km^{3}/day, most of which occurred within the first 7 days, contributing to a total water loss of ~58 ± 6 km^{3} (1σ; equivalent to 61% of the total deposited stormwater amount). We can also obtain a first-order estimate of groundwater flow by separating the *SG* component using an empirical hydrograph separation filtering technique, which assumes that groundwater flow is a low-frequency pulse (*24*). This indicates that groundwater flow contributed to ~12 km^{3} (14%) of overall water loss, at ~1.25 km^{3}/day (red line in Fig. 4A). Comparing the *SG* component with the measured equivalent from USGS streamflow data (cyan versus orange lines, respectively; Fig. 4A) shows agreement in basin lag times (of 1 to 2 days, that is, time lag between centroid of rain and centroid of *SG*) but that *SG* has almost four times greater total volumetric water loss in the first 7 days following Harvey’s first arrival. This suggests that this initial, large pulse of stormwater loss likely reflects the “ungauged” component, not captured from the available 39 river gauges (fig. S13C for location of river gauges), and most likely flowed downstream of available river gauges, and/or drained directly into the Galveston Bay adjacent to Houston and the ocean. Our results indicate that GPS estimates of water storage can provide observational constraints on the characteristics and response of a region’s hydrologic system to extreme, transient precipitation events, from its storage capacity, the initial rapid loss of stormwater that is difficult to comprehensively capture via traditional streamflow gauges, and to the gradual recovery time of the system.

The amount of standing surface water (a component of TWS) that later reaches nearby streams via surface runoff poses a continued flood hazard to downstream rivers and dams (*3*). Short-term operational flood forecasts use hydrologic and hydraulic models to estimate flooding potential and changes of river levels over 1- to 2-day time scales during and following rainfall events (*25*). Some hydrologic models have assimilated remote-sensing data sets that can measure the degree of preflood soil saturation, which helps constrain the amount of rainfall that infiltrates versus that which leads to surface runoff (*26*, *27*). These satellite measurements include monthly TWS estimates from GRACE gravimetry, which constrains the degree of basin saturation (*26*), and soil moisture content from MetOp (Meteorological Operational satellite programme) and SMOS (Soil Moisture and Ocean Salinity) (*27*, *28*), both of which have been shown to improve the predictive skill score of medium- to long-range forecasts (day-month time scales) of runoff and streamflow levels (*26*–*28*). Combining preflood surface saturation measurements from satellites with daily water storage information from GPS during flooding can provide additional constraints for hydrologic models on the evolving state of deposited water and its probability of runoff (versus infiltration), potentially improving the reliability of short-term forecasts of river levels during and following large rain events. However, for GPS to be used in this manner would require near–real-time positioning with 1-day latency. To assess the feasibility of this approach, we reprocessed the GPS positions using satellite orbital and clock information that would be available at 1-day latency (producing daily static positions; see Methods for details and fig. S14). We found that comparison of the rapidly determined GPS positions with those derived from final precise orbits and clocks (available after 7 days; Figs. 1 and 2), have equivalent precision (WRMS repeatability of the unfiltered time series of 6.16 and 6.35 mm for the rapid and final time series, respectively) and strong consistency (median ρ = 0.85). In addition, we find that the ICA filter method is also able to extract the hydrologic signal from the rapidly available data (fig. S15). Although near–real-time analysis of GPS for water monitoring and forecasting will be challenging to implement, requiring removal of systematic errors and sufficient coverage of continuous GPS stations, TWS estimated from the final precise data (available after 7 days and presented here; Fig. 2) could have practical use in helping calibrate or validate hydrologic models that simulate previous rainfall-runoff events for helping mitigate the effects of future floods. Daily estimates of water storage from continuous GPS networks show potential for use in near–real-time flood alerts, monitoring and forecasting for better understanding how downstream river and dam levels could change in response to direct rainfall and stored water in the days during and following major storm events, information that would be vital to flood managers (*29*, *30*).

Our analysis of GPS surface motions during the passage of Hurricane Harvey demonstrates—for the first time—that it is possible to robustly quantify daily changes in transient water storage following large hydrologic events. This underscores a potentially useful role that dense and continuous GPS networks could play in monitoring flooding risks and guiding the necessary response and mitigation. The global availability of continuous GPS networks provides the potential for near–real-time, remote monitoring of water storage extent in the days during and following large hurricanes and other types of storms, where TWS is not only a key component of the hydrologic cycle but also a potential, continued flood hazard once the storm system has passed. As the Gulf of Mexico experiences warmer air and water temperatures in an ongoing warming climate, the occurrence of large influxes of water from powerful hurricanes like Harvey is unlikely to decrease (*31*, *32*). These climatic conditions therefore warrant addressing the efficacy of current drainage networks and the development of innovative flood mitigation strategies in rapidly expanding urban areas such as Houston, where the type of analysis presented here can aid in preparing and monitoring for the effects of flooding in future similar events.

## METHODS

### Processing of GPS data

Daily changes in TWS were estimated using the vertical and horizontal components of the GPS time series that recorded surface deformation due to changes in hydrological loading. We estimated the daily GPS positions using Jet Propulsion Laboratory (JPL) GIPSY-OASIS II in precise point positioning (PPP) mode (*33*), with resulting position solutions aligned to the ITRF08 reference frame. To correct for tropospheric delays, we used the Vienna Mapping Function 1 (VMF1) and nominals that are estimated every 5 min (*16*, *17*), and we made corrections for ocean tidal loading (*34*). In postprocessing, we corrected for nontidal atmospheric and nontidal ocean loading using corrections from (*10*), which included hurricane-related pressure changes and the storm surge. The reduction in variance from the nontidal ocean and atmospheric corrections range up to 10 and 21%, respectively. Before these corrections, we applied a criterion to remove stations that exhibited large noise in their time series or were missing extensive data. We removed 19 stations that were missing more than 30% of their positions. We then removed an additional 15 stations that exhibited a variance that was three times that of the median variance estimated from the entire network. In the case that a particular GPS time series was missing a segment of data (those that passed the 30% criterion), we followed the approach of (*13*), where <3 consecutive days of missing data in the time series were replaced using a linear interpolation, and for data gaps of ≥3 days, we replaced missing values using the average position of the entire network. We found that the number of missing data positions on a given day for all stations across the entire network was <2% (equivalent to <4 of 219 stations; see fig. S16). This indicates that, when taking the network average position (in the occurrence a station was consecutively missing more than 3 days of data), it was not obtained from a small subset of the GPS network, and is therefore unlikely to be erroneously mapped into CME and removed by the ICA filter.

### Spatiotemporal filtering—ICA

Principal component analysis (PCA) and ICA are forms of blind source separation techniques that seek to unmix data to identify the underlying and unknown sources (*13*). Here, the goal of using these feature extraction techniques was to isolate and remove CME that appeared as spatially correlated noise from the GPS data and to separate it from the hydrologic load signal. PCA is a statistical approach that uses an orthogonal linear transformation to reproject the data from a set of possibly correlated variables to a set of linear uncorrelated variables of maximum variance, called principal components (PCs). The PCs, which represent the temporal basis functions of the data, were generated by projecting the data onto a set of orthonormal basis vectors that were derived from an eigenvalue decomposition of the data covariance matrix. These eigenvectors describe the spatial pattern of the data and are referred to here as spatial responses. These are ordered by the percentage of variance explained, where the first component denotes the source of highest variance and contributes the most motion to the GPS network. Because the covariance matrix of the data matrix (*n* GPS stations by *t* time samples) is full rank, the eigendecomposition provides *n* set of PCs and *n* set of spatial responses. However, because most GPS time series do not follow a Gaussian distribution (that is, the underlying process is not Gaussian) (*15*) and PCA uses second-order statistics (variance) and assumes that the underlying sources are Gaussian, it is not an optimal method and is susceptible to mixing sources across different components.

ICA is similar in concept to PCA; however, it finds sources of maximum independence instead of minimum correlation—the former a stronger condition as not all functions that are uncorrelated are independent. ICA is advantageous over PCA as it uses higher-order statistics (fourth order; for example, negentropy in this case) and assumes that the underlying components are non-Gaussian and statistically independent.

To perform the ICA analysis, we organized the data matrix (*X*_{n×t}) into *n* rows and *t* columns, with each element representing displacement in a certain direction (for example, vertical component of GPS) measured by the *n*th GPS station. ICA was applied separately to the vertical, east, and north components. For each row, we subtracted the sample mean and then whitened the data. As shown in Eq. 1, the observed data matrix was assumed to be some transformation (*Q*_{n×r}) from a set of *r* unknown time-varying sources (*S*_{r×t})(1)

The task of ICA is to determine the unmixing matrix (*W* = *Q*^{−1}) so that we could determine the underlying the sources (*S*) from the data (*X*). In PCA, the *Q* matrix is a linear orthogonal transformation that maximizes the variance of the rows of *S*, while in ICA it uses a similar generalized form of Eq. 1, but *Q* is a nonlinear transformation that maximizes the statistical independence of the rows of *S*. Here, we used the reconstruction ICA (rICA) approach (*9*) to estimate the unknown sources. This differs from other ICA methods such as fastICA (*35*) by swapping the orthonormality constraint applied to *W* (that is, *WW*^{T} = *I*), with a reconstruction penalty term added explicitly to the objective function, giving the benefit of using unconstrained solvers [see equation 2 of (*9*)].

### Number of appropriate ICs

One limitation to ICA approaches is that they cannot determine the number of appropriate components to decompose the data or their ordering. The number of components chosen is important as one too many may result in incorporating noisy sources and too few may mix the signals together. To determine the appropriate number of components, we used two different stopping methods applied to the eigenspectrum derived from PCA (*13*).The first stopping method is a measure of separability of the eigenvalues (λ_{i}), where we used a standard rule of thumb (*36*) to assess which eigenvalues (λ_{i}) exceed that expected from a random process. If the separation between eigenvalues (Δλ_{i} = λ_{i} − λ_{i−1}, for *i* > 1) falls below the uncertainty ∂λ_{i} = λ_{i}(2/*n*)^{2}, then the component becomes more difficult to separate from its neighbor and from noise. The second approach uses Horn’s parallel analysis (*37*), a Monte Carlo simulation approach that randomly scrambles the data to create a suite of random samples and its simulated eigenspectrum and their uncertainties. If the observed λ_{i} exceeds the 95% confidence interval of the simulated λ_{i}, then the component is retained. We found that both methods indicate that four PCs are significant, and we used this number to decompose the data.

### Ordering of ICs

Because of the whitening of data as a preprocessing step in rICA, the ordering of the components cannot be determined directly. To resolve this, we estimated the ordering using a ratio following (*38*), which characterized the contribution of the *i*th IC to the observed unfiltered motions (*X*)(2)

The importance (*H*) of each IC (*i* = 1,…,4) was estimated at each station (*k*), as the ratio of the observed motions (*X*) with time (*t*), with the predicted motion due to the *i*th IC at that station. For each IC, we then estimated the median of the ratios across all stations to find the relative contribution that each component made to the overall network, with lower values indicating higher importance. As discussed in text S1, the IC we associated with CME was found to have the lowest median *H* value (that is, the most important), consistent with the PCA that found that the same CME component had the largest eigenvalue.

### Identifying CME

To determine which component corresponds to CME, we followed the criteria of (*13*). A component is considered to represent CME if >50% of stations exhibit a significant normalized spatial response (>25% of the maximum) and the eigenvalue exceeds 1% of total variance. In addition, we inspected the temporal pattern of each component to assess whether it exhibits changes corresponding to the timing of Hurricane Harvey, and estimated its spectrum, where previous studies have found CME to follow a flicker noise process (slope of −1) (*15*).

To determine CME from the GPS data, we initially applied ICA to the entire network (fig. S17). However, we found that this absorbed some of the real hydrologic signal (fig. S2, A and B). Instead, we estimated CME from a set of reference stations (*n* = 21) located in the northwest of the network, which are close enough to stations affected by Harvey but distal from regions of known precipitation (fig. S1) (*20*). Therefore, this approach is equivalent to most geodetic studies that take GPS positions relative to a group of reference stations distal to some target source (*11*).

CME was identified on IC1 (fig. S1) following the criteria described above, and which we found followed flicker noise (fig. S2C), consistent with a previous study of 259 GPS data in China (*15*). After this, we removed it from the rest of the network by simply multiplying the median spatial response of the 21 reference stations with IC1 and subtracting it from the time series of each GPS station. We found that, when applying ICA again, there was no remnant “CME component” found, indicating that it had been successfully removed (fig. S3).

For the horizontal component of the GPS motion, we followed the same approach as the vertical component, where we estimated CME from the same subset of stations and removed this from the data. We then reapplied the ICA analysis to extract the hydrologic loading signal. Similar to the vertical component of motion, we found two ICs that exhibited spatiotemporal changes consistent with the timing and spatial pattern of Harvey. For both the east and north components of the GPS motion, we interpreted the third and fourth ICs as reflecting hydrological loading (figs. S4 and S5), with clear deviations in time that corresponded to both landfalls of Harvey and spatial responses centered around Houston and west Louisiana.

### Inversion of GPS for water loading

TWS was resolved as daily changes in water mass on a grid of 0.25° × 0.25° nodes. Changes in the GPS-derived vertical and horizontal components were mapped to a surface mass load using Green’s functions for a spherical self-gravitating layered Earth model (*4*, *19*). We minimized the objective function in Eq. 3 for the vertical and horizontal elastic responses to the daily water mass changes, with spatial and temporal regularization terms, using a positively bounded L_{1} norm inversion following (*39*)(3)where *G*_{w} and *d*_{w} are the weighted Green function matrix and data vector, respectively, and other terms are defined below. *G*_{w}*m*_{t} = *d*_{w} together with *Sm*_{t} = 0, and *Um*_{t} = *Um*_{t−1} can be expanded as follows(4)where *W* is a diagonal weighting matrix associated with the formal daily uncertainty of the GPS position and *G* are Green’s functions (*n* × *p* matrix) relating *p* number of water mass patches loading a spherically layered elastic Earth (*19*) to surface deformation observed at *n* number GPS stations, with subscript *v* and *u* denoting vertical and horizontal components, respectively. is a vector of length *n* containing the vertical component (*v*) of the observed GPS displacement (*d*) on day (*t*), *m*_{t} is the water mass that we solve for each day (*t*), and superscripts *e* and *s* denote east-west and north-south component, respectively. The solution was also regularized with temporal and spatial smoothing to suppress large and unrealistic variations of water mass between patches in space and time. λ and β are time-invariant coefficients that control the strength of spatial and temporal smoothing, respectively. *S* and *U* are the smoothing operators that estimate the current model’s (*m*_{t}) spatial curvature and the difference from the previous model (*m*_{t−1}), respectively, achieved with a central and backward finite-difference approximation, respectively. The choice of these smoothing factor values is discussed later below. The above problem was solved for using an outlier-resistant L_{1} iterative solver following the method of (*39*)(5)where *M* is the leftmost term in Eq. 4, *M*^{T} is the transpose (which contains the weighted *G* and the regularization terms), and *b* is the rightmost term in Eq. 4 (which contains the weighted data vector, with the regularization terms appended below). *R* is a diagonal weighting matrix, with diagonal elements that are the absolute values of the reciprocals of the residuals(6)(7)Therefore, *R* is a nonlinear function of *m*_{t}, and the nonlinear system in Eq. 5 was solved for using an iteratively weighted least squares solver algorithm (*19*) that repeats the inversion for a new *R* once the model and residual vectors converge. The initial solution used to estimate *R* was taken from a positively bounded L_{2} norm inversion *m*_{t} *=* (*M*^{T}*M*)^{−1}*M*^{T}*b*.

The horizontal data were equally weighted in the inversion as the vertical data. Although the data uncertainty is smaller for the horizontals than for the verticals (on average 2 mm for horizontal and 5 mm for vertical at 1σ), meaning slightly higher weights in *W*, water mass changes are primarily sensitive to vertical motions. Because mass loading primarily causes vertical surface motion, the Green’s function values that relate the two are almost an order of magnitude larger, and therefore, water mass changes are more sensitive to vertical GPS changes than to the horizontal GPS changes. This explains the minor difference in the solution when including the horizontal in the inversion. However, we noted that the main effect of including the horizontal data is that it appears to improve the fidelity of the solution, concentrating water mass and giving a less smoothed result. This is not surprising, given that the horizontals are sensitive to local mass changes. Uncertainties of the TWS solution were obtained from 10,000 Monte Carlo simulations of the GPS data given the daily uncertainty (fig. S18).

### Validation of model and filtered GPS data

To attempt to validate our GPS data, we compared our filtered time series to surface motions predicted from a hydrologic model. For each day, we forward model the daily averaged loading on a spherically layered elastic crust from TWS simulated by the NLDAS (fig. S19) (*18*). NLDAS estimates hourly changes in TWS at 1/8° spatial resolution over the continental United States, driven by observed precipitation. We estimated mean daily TWS from the NLDAS model as the change relative to the mean simulated TWS from 8 days before Harvey, therefore isolating the effects of water variation from Harvey itself.

Soil moisture measurements from NASA’s Soil Moisture Active Passive (SMAP) satellite (fig. S20) also show qualitatively good agreement in distribution of the saturated surface with our estimates of TWS (Fig. 3). We note that large soil moisture and TWS as simulated from the NLDAS model were also found northwest of Houston, as seen in our inversion result (Fig. 3) and highlighted by the red circle in fig. S18.

### Stream discharge measurements

Subdaily river discharge measurements were obtained from 39 rivers along the Gulf Coast from the USGS (https://waterdata.usgs.gov/nwis/uv), using approved data. For each river, we selected discharge data from the closest available station to the mouth of the river to maximize the catchment area (fig. S13). The total water volume loss due to river discharge with time, which is shown in Fig. 4, was estimated by summing the integrated discharge profiles with time for all rivers and removing a base flow determined as the average discharge before Harvey. We note that this removal of the average discharge is not an attempt to remove baseflow (which can take on various methods) but rather to capture the deviation of total river discharge (and volume) from its typical flow. Rivers that are missing portions of data in the time series were corrected for by using a shape-preserving, piecewise cubic Hermite polynomial interpolation. Although we have obtained measurements of discharge from all the major rivers (8 major rivers) and 31 smaller ones, we are likely missing smaller outlets and consider the estimated river discharge value a minimum estimate.

### Poroelastic effects

To investigate whether poroelastic effects from recharge of groundwater bias the GPS motions, we analyzed 30 time series of wellhead levels obtained from the USGS around the Houston area and compared these to nearby GPS stations. Recharge of aquifers leads to elastic expansion and uplift that could be detected at nearby GPS stations. Here, we analyzed wellhead levels to estimate how much water diffused into the aquifer and increased the water level during the period of observed GPS uplift. If poroelastic expansion does affect GPS motions, then we would expect largest water increase during largest GPS vertical uplift rates. GPS uplift occurred over a period of ~5 weeks, however, with highest rates within the first 4 days.

We found that the comparison of most groundwater levels with vertical GPS displacement are inconsistent with surface motions expected from poroelastic effects. Almost all groundwater sites (87%) either show no deviation in behavior during Harvey or exhibit water level increases during largest GPS subsidence, which is followed almost immediately by a groundwater drop, when GPS uplift is observed. This surface motion is the opposite of that expected from poroelastic deformation. From the 30 groundwater levels, we only observed 4 that show an increase in water level during GPS uplift. However, we note that changes of water levels before and well after Harvey, which are similar in amplitude and duration (days-weeks) to water changes during Harvey, showed little or no response from the GPS stations. This would suggest that, from the current data, GPS stations are not sensitive to groundwater level changes during recharge over day-week time scales.

We estimated vertical displacement expected from poroelastic aquifer recharge using the wellhead data and assuming an elastic clay storativity that represents the Chicot and Evangeline aquifers (*40*). These estimates were then compared to the nearest GPS stations that were no further than 25 km away (fig. S6). To estimate vertical displacement of the aquifers (Δ*d*) due to changes in wellhead level (Δ*h*), we assumed a storativity coefficient (*T*) of 2 ×10^{−3} derived from (*40*), where Δ*d* = Δ*h* × *T*. From this more direct comparison, we found that the observed GPS motions were inconsistent with that expected from poroelastic motion from aquifer recharge (fig. S6). We found that, in almost all cases, the sign, amplitude, and functional form disagrees between the observed surface motions and that expected from poroelastic aquifer recharge. Instead, the aquifer levels can be interpreted to reflect direct changes in water across at the site and validate the filtered GPS time series.

### GPS processing for near–real-time analysis

We reprocessed the GPS time series of the 219 stations to assess whether the position precision would be sufficient to resolve the hydrologic loading. We generated a time series of the daily static position using GIPSY-OASIS II in the PPP mode, with JPL rapid orbits and 5-min clocks. These are available daily by 12:00 UT, allowing processing of data from the previous 24-hour UTC day. PPP configuration is the same as for the final time series but without second-order ionosphere corrections, as is recommended because rapid orbits are generated without these. In addition, the GPT2 troposphere model was used for mapping functions and nominal values, because VMF1 information is not available in time for rapid processing.

The comparison of the GPS rapid solutions versus the precise (the latter used in this study) shows strong agreement (fig. S14), with a median ρ of 0.86 for all stations and a median percent of variance reduction (eq. S1) of 75%. The degree of precision, or daily repeatability (estimated using a WRMS), is also similar between the two, with medians of 6.16 mm for rapid solutions and 6.36 mm for precise. Furthermore, we then applied the same ICA method as that used for the precise solutions, finding a similar result in both the temporal and spatial behavior (fig. S15), with both landfalls separated onto different components. We note that the ICA method does not require the entire time series, and fig. S15 shows ICA applied to the rapid time series after the first day of initial landfall of Harvey, indicating that it can successfully capture the hydrologic signal.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/4/9/eaau2477/DC1

Materials and Methods

Text S1. Analysis of CME

Text S2. Smoothing parameters

Text S3. Synthetic test

Fig. S1. Estimating CME from subset of stations.

Fig. S2. Comparison of CME estimate from two regions of GPS network illustrating difference in degree of mixing of the hydrological signal.

Fig. S3. ICA analysis of vertical positions after CME removal.

Fig. S4. ICA analysis of horizontal positions (north component) after CME removal.

Fig. S5. ICA analysis of horizontal positions (east component) after CME removal.

Fig. S6. Comparison of vertical motions from GPS and expected poroelastic rebound.

Fig. S7. Comparison of vertical motions from GPS and predicted from hydrologic (NLDAS) model.

Fig. S8. Comparison of north motions from GPS and predicted from hydrologic (NLDAS) model.

Fig. S9. Comparison of east motions from GPS and predicted from hydrologic (NLDAS) model.

Fig. S10. Percent of variance reduction (POVR) for different temporal and spatial smoothing values.

Fig. S11. Accumulation and dissipation of a synthetic water volume source as a function of time in days.

Fig. S12. Inversion result of synthetic test.

Fig. S13. River discharge and cumulative water profiles.

Fig. S14. GPS vertical positions processed using rapid orbits to assess use for near–real-time applications.

Fig. S15. ICA filtering of GPS solutions processed using rapid orbits, which were available with 1-day latency.

Fig. S16. Raster matrix highlighting space-time distribution of missing data points for all GPS time series.

Fig. S17. Components of spatial and temporal variations of GPS data separated using ICA.

Fig. S18. Model uncertainty (1σ) of each grid node in water volume (km^{3}).

Fig. S19. Daily simulated TWS from NLDAS hydrologic model.

Fig. S20. Soil moisture measurements from NASA’s SMAP satellite.

Fig. S21. TWS inversion estimates for rest of days not shown in Fig. 3.

Data file S1. GPS raw and filtered time series positions shown in Figs. 1 and 2.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**We would like to thank the three anonymous reviewers for their helpful suggestions. We also thank J. Fisher for providing the ET data and J. Reager, T. Farr, and A. Gualandi for helpful discussions.

**Funding:**Part of this research was supported by the NASA Earth Surface and Interior focus area (ESI) and performed at the JPL, California Institute of Technology. Funding for this project was provided under a NASA Postdoctoral Program fellowship to C.M. administered by the Universities Space and Research Association through a contract with NASA, and NASA ESI grant NNX17AE01G awarded to R.B. and Y.F.

**Author contributions:**R.B. contributed to formulation of scientific analysis. C.M. performed postprocessing of GPS data and developed the numerical inverse scheme. K.M. helped acquire supporting data. A.W.M. helped with processing of GPS time series. Y.F. and S.A. helped with inversion modeling for water loading. D.B. helped with inversion scheme. D.F.A. helped supply supporting data. All authors participated in manuscript revision.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**Precipitation data can be accessed from www.emc.ncep.noaa.gov/mmb/ylin/pcpanl/stage4/. USGS river discharge can be accessed from https://waterdata.usgs.gov/nwis/uv. This material is based on data provided by UNAVCO through the Geodesy Advancing Geosciences and EarthScope Facility with support from the NSF and NASA under NSF Cooperative Agreement No. EAR-1261833. GPS raw and filtered time series are available as auxiliary data and can be found at a data repository from https://chrismillineractivetectonics.com/data-codes/. Barker and Addicks reservoir data can be found from (https://waterdata.usgs.gov/nwis/inventory/). All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors. We would also like to thank the Louisiana spatial reference center (LSRC) at the Louisiana center for geoinformatics (C4G) for access to their data.

- Copyright © 2018 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).