A total of 18 Corpus Christi Bay surface water samples were collected within 10 m of the waterline at University Beach beginning on 23 August 2017 and continued until 4 February 2018. Henceforth, we refer to surface water samples from Corpus Christi Bay as surface water samples. All surface water samples were in line with the groundwater and porewater samples ( Figure 1C ). Two pre-hurricane surface water samples were collected 2.5 and 1.7 days before the first landfall. Post-hurricane, 11 samples were collected daily beginning on 30 August, the fourth day after the mainland (and second) landfall, and three more samples were collected weekly between 16 September and 30 September. An additional two surface water samples were collected in February 2018. Two porewater samples were collected pre-hurricane landfall and one in February 2018. In addition, six groundwater samples were collected from two nested monitoring wells near the monitoring site ( Figure 1C ): one developed between 5 and 6 m in the unconfined aquifer (shallow) and another between 10 and 12 m in the semiconfined water producing unit (deep). Groundwater samples were collected from each formation 2 days pre-hurricane, post-hurricane on 20 September, and in February 2018. All samples were analyzed for inorganic nutrients, dissolved organic carbon (DOC), total dissolved nitrogen (TDN), and DOM molecular characterization. All samples were collected using standard techniques for surface water, porewater, and groundwater (Douglas et al., 2021a).

Nitrate (NO₃ ⁻), nitrite (NO₂ ⁻), and ammonium (NH₄ ⁺) were measured on a Seal QuAAtro nutrients autoanalyzer (EPA, 2011) at the Geochemical and Environmental Research Group (GERG), College of Geosciences at Texas A&M University College Station, with method detection limits (MDLs; in µmol L⁻¹) of 0.11, 0.012, and 0.057, respectively. The MDL is defined as Student’s t for 99% confidence level times the standard deviation of seven replicate measurements of the same low-level sample or spiked sample. Measurements of DOC and TDN were conducted using a Shimadzu TOC-V analyzer with a nitrogen module, with an MDL of approximately 1 mg·L⁻¹. Dissolved organic nitrogen (DON) was estimated as the difference between TDN and dissolved inorganic nitrogen (DIN = NO₂ ⁻ + NO₃ ⁻ + NH₄ ⁺).

Samples for DOM molecular characterization were collected, processed, and analyzed with UPLC-OT-FT MS as described in detail in Douglas et al. (2021a). Additional OT-FT MS methods and specific OT-FT MS instrument settings can be found in the Supplementary Material . Briefly, samples were collected in 1-L polycarbonate amber bottles previously cleaned with Alconox and HCl solution and rinsed with the sample. Immediately after collection, samples were sterile filtered through 0.2-μm polyethersulfone filters conditioned with 5–10 ml of sample to remove particulates and frozen until processed. Prior to extraction, 100 ml of each sample was acidified to pH 2 with HCl to increase extraction efficiency for organic acids and phenols. DOM samples were then prepared for mass spectrometric analysis using 200 mg, 3 ml Bond Elut-PPL solid phase extraction cartridges following the procedure recommended by Dittmar et al. (2008). Preceding elution of DOM, two cartridge volumes of 0.01 HCl were passed through the cartridge to completely remove salts and then were dried. The samples were eluted with 6 ml of Optima LC/MS-grade methanol. The methanol eluents were dried in a Centrivap benchtop concentrator and reconstituted with 1 ml of 95:5 water to acetonitrile. 

Aliquots of the sample were analyzed on a 1.7-μm ACQUITY UPLC BEH C18 reversed-phase column (Waters, 30 Å, 1.7 µm, 2.1 mm × 100 mm) via a heated electrospray ionization (H-ESI) source to Orbitrap Fusion Tribrid Mass Spectrometer (Thermo Scientific, Waltham, MA, USA) operated under positive mode. The positive mode was used for its preferential ionization of nitrogen-containing organic matter to facilitate comparison with and discussion of changes in inorganic nitrogen species. Between each sample analysis, an 8-min column re-equilibration was performed. For MS², the isolation window was set at 0.7 m/z while performing both collision-induced dissociation (CID) and higher-energy collisional dissociation (HCD) using an ion trap mass spectrometer as the detector.

Compound Discoverer software 3.1 (Thermo Fisher) was used to identify the DOM compounds according to the following conservative criteria: 1) a signal-to-noise ratio of above 3, 2) a minimum of five mass scans per chromatographic peak, 3) a minimum peak intensity of 50,000, and 4) at least one isotope peak detected. All chromatography spectra retention times (RTs) were aligned using an adaptive curve with a maximum shift of 0.2 min and 5-ppm mass tolerance. The ratio of the M+1 to parent peaks was used to confirm the number of carbon atoms, and the M+2 peak ratio was used to confirm the presence of the S atom in the compound. If multiple adducts were detected ( Supplementary Material ), all adducts of the same compound were grouped together with an RT tolerance of 0.2 min. Based on our observed analysis, there were limited cases where multiple adducts were detected for the same compounds. A total of 1,427 out of 1,822 compounds detected had a single adduct ([M+H]⁺¹), and the maximum number of adducts detected for a single compound was 6.

A molecular formula calculator (Molecular Formula Calc version 1.0 NHMFL, 1998) was used to determine empirical formula matches with the following parameters: C_4–100H_4–200O_0–50N_0–10S_0–3P_0–3. Molecular formulas that are unlikely to occur in nature or are not chemically possible were removed, as described in detail in Abdulla et al. (2013). Internal standards were not used when analyzing these samples; thus, the calculated masses of the assigned formulas are all within 2.0 ppm of the masses detected by OT-FT MS.  For a more detailed description of these analyses, the reader is directed to Douglas et al. (2021a).

The data were gap filled and normalized to constant sum in Compound Discoverer and exported for statistical analyses. Compound presence or absence in each sample was determined if it met or did not meet the four criteria listed above. For absent compounds, a gap filling approach that integrated the area where the peak was expected but not detected was performed to prevent skewing toward missing values (i.e., zero values) in the statistical analyses. No bias is introduced by integrating these gap-filled areas, as they usually correspond to spectra noise. Peak magnitude-weighted O/C_w, H/C_w, and N/C_w were calculated, as described in Douglas et al. (2021a), to characterize each sample and facilitate discussion of DOM composition (Sleighter et al., 2010; Schmidt et al., 2017).

A molecular lability boundary (MLB), developed from FTICR-MS molecular data visualized in van Krevelen diagrams with chemical class regions, separates the data into more and less labile constituents. This is based on the higher degrees of hydrogen saturation of the protein-, amino sugar-, and lipid-like region linkage to more bioavailable fractions of carbon (D'Andrilli et al., 2015). Constituents of DOM above the MLB at H/C ≥ 1.5 correspond to more labile material, whereas constituents of DOM below the MLB (H/C < 1.5) exhibit less labile, more refractory character (grouping in more lignin-, tannin-, and condensed aromatic-like regions of the van Krevelen diagram). The percentage of more labile constituents (MLB_L) was determined by the number of molecular formulas with H/C ≥ 1.5 over all O/C (0.0–1.2) divided by the total number of molecular formulas in the sample or by the number of formulas in a heteroatom group, multiplied by 100. Classification of less labile, or more refractory components (MLB_R), from each carbon source was calculated similarly for molecular formulas below the MLB (H/C < 1.5). These calculations are defined as chemical species richness for labile and refractory contributions. The diversity of chemical species was calculated by incorporating the normalized peak height of each mass spectral peak to determine the weighted effect of each ionized elemental composition on the overall distribution in each group (i.e., labile: MLB_wL or refractory: MLB_wR). Although chemical class regions in van Krevelen diagrams contain a degree of ambiguity, and the degree of lability can also be a function of oxygenation; this conservative MLB calculation can provide an unbiased way to distinguish between more or less labile OM by high-resolution MS (D'Andrilli et al., 2015).

Cluster analyses and principal component analyses (PCAs) were performed to assess the temporal variability of DOM composition using SAS 9.4 software. Cluster analyses were performed using PROC CLUSTER with Ward’s minimum-variance method after standardizing DOM peak magnitudes with PROC ACELUS, which transforms the data such that the resulting within-cluster covariance matrix is spherical when you suspect the data contain elliptical clusters without known cluster membership or number of clusters, and PROC STANDARD, which creates z-scores for each variable to a mean of 0 and standard deviation of 1. The number of clusters was determined from a visual inspection of the dendrogram and the values of the cubic clustering criteria, pseudo-F, and t² statistics. To explore how sample groups differed based on compound class and structure, the O/C and H/C ratios of the assigned molecular formulas for compounds present in each cluster were plotted in van Krevelen diagrams with generalized compound class regions identified (Koch and Dittmar, 2006). The differences in DOM composition between the clusters aid in identifying the potential sources and transformations of compounds in the time series as well as potential environmental drivers.  

PCA and van Krevelen diagrams were used to visualize the correlation between multivariate data in the samples. For PCA, the data were standardized to a normal distribution using PROC STANDARD so that scales were the same for all variables. PCA was performed using the PROC FACTOR procedure in principal components mode with VARIMAX rotation for axis rotation (Douglas et al., 2021a). A sample score plot and a variable loading plot, which shows the distribution of the samples and variables (i.e., individual compounds), respectively, were generated.

Spearman’s rho coefficient (r) was used for correlation analyses of nutrients and hydrologic parameters (i.e., water level and creek discharge). This type of analysis is best suited for parameters not normally distributed and are not expected to have a linear relationship, thus making Pearson’s correlation coefficient unreliable (Hauke and Kossowski, 2011). For this study, the 95% confidence interval (or p-values ≤0.05) is deemed as an indication of a significant correlation.

Concentrations of NO₃ ⁻, NO₂ ⁻, and NH₄ ⁺ were substantially greater following landfall and the accompanying ~0.5 m drop in water levels from before the storm ( Figure 3B ). However, only surface water NH₄ ⁺ had a strong negative correlation with water level (p-value = 0.03), while NH₄ ⁺ and NO₂ ⁻ were positively correlated (p-value = 0.02 and p-value < 0.01, respectively) and DOC was negatively correlated (p-value = 0.04) with Oso Creek discharge ( Table 1 ). In contrast, surface water DON concentrations were slightly higher before the hurricane landfall and initially decreased post-landfall ( Figure 3C ). Concentrations of DON, NO₃ ⁻, NO₂ ⁻, and NH₄ ⁺ fluctuated with a decreasing trend corresponding with the most negative surge, with water levels decreasing to MLLW or lower multiple times over the 4 days after the second landfall. The most negative surge followed 2 days after Harvey’s final landfall, though it occurred approximately 280 miles (451 km) northeast of Corpus Christi in southwestern Louisiana. Though concentrations of DON and NH₄ ⁺ initially followed the same trend, they began to fluctuate with inverse trends to each other starting on day 9 while returning to MSL ( Figures 3B, C ). As DON spiked with a slight drop in salinity in the early days post-landfall, NH₄ ⁺ decreased ( Figures 3A, C ). Concentrations of NO₃ ⁻ remained relatively low but peaked along with DON slightly on day 9 and substantially on day 13 and were overall positively correlated to each other (p-value < 0.01). Concentrations of NH₄ ⁺ peaked slightly on day 10 with NO₂ ⁻ concentrations peaking over the following 2 days before the substantial NO₃ ⁻ peak on day 13, indicating a recent local pulse of reduced inorganic N followed by nitrification processes.

All nutrients and DOM concentrations stabilized on the lower end for the following 3 weeks, with some slight differences. With a steady overall increasing trend in water levels following day 5, DON and NH₄ ⁺ seem to have been modulated to a certain extent by salinity fluctuations (p-value < 0.01) and the daily tidal cycle (p-value = 0.03), respectively. Concentrations of DOC were overall lower post-landfalls. Two periods of decreasing concentrations correspond to declines in salinity, which may be pulses of freshwater inflows.

The groundwater and porewater signatures of different dissolved inorganic and organic nutrients were used in this study, representative of pre- and post-landfall conditions, as means to evaluate the potential for groundwater–surface water interaction. Groundwater in the nearby island water table (e.g., unconfined aquifer) was characterized by low to brackish salinity (0.37–6.47, x̄ = 3.1), both pre- and post-landfall, while for the semiconfined (deeper) aquifer, salinities were generally hypersaline (30–54, x̄ = 44.3) though slightly fresher post-landfall ( Figures 1C , 4 ). Overall, the salinity of shallow groundwater was lower than that of porewater, surface water, and deep groundwater, likely due to recharge from local rainfall, irrigation, and limited connection with the estuary. Similarly, the signatures of different dissolved inorganic and organic nutrients reveal shallow groundwater as chemically distinct from deeper groundwater ( Figure 4 ).

While pre- and post-landfall changes in pore- and surface water salinity had trends like the deep and shallow groundwater (e.g., a shift toward lower salinity following the landfall), changes in NO₃ ⁻, NO₂ ⁻, and NH₄ ⁺ concentrations were not consistent ( Figures 4A–C ). The concentrations of NO₂ ⁻ and NO₃ ⁻ were the highest in shallow groundwater a few months after the landfall, when the system was more hydrologically stable. Deeper groundwater had a range of NO₂ ⁻ concentrations that encompass porewater and most of the surface waters ( Figure 4B ). Concentrations of NH₄ ⁺ were the highest in porewater a few months post-landfall and about threefold lower pre-landfall, with levels approaching pre-landfall shallow groundwater ( Figure 4C ). However, NH₄ ⁺ in porewater exceeded deeper groundwater and post-landfall shallow groundwater. A trend with salinity for groundwater was not obvious for inorganic N, except for NO₃ ⁻, which was the highest at the low salinity end, in shallow groundwater ( Figure 4A ).

Salinity and DON ( Figure 4E ) exhibit a negative trend with surface water and porewater falling between the shallow and deeper groundwater fields. However, salinity and DOC ( Figure 4D ) showed a positive trend in each source, with the pre-landfall samples at the high end and samples from months after landfall at the low end, which may be due to a seasonal effect, but the sample size is too small to speculate further. The largest DON concentrations were in the shallow groundwater and the lowest in the pre-landfall deeper groundwater. Deeper groundwater DON was higher post-landfall, while the opposite was observed for the shallow aquifer. In comparison, pre-landfall porewaters were more enriched in DON than surface water and deeper groundwater. These trends in DON among the different environments and pre- and post-conditions reveal a mixing pattern between the two aquifers and porewater, and likely exchange with surface water.

A total of 1,821 unique compounds were identified by DOM molecular analyses on pre- and post-hurricane surface water and groundwater samples (n = 26). Of the total unique compounds, 1,452 (80%) were assigned formulas that were primarily CHO (37.9%), CHON (39.6%), CHONS (11.1%), and CHOS (4.9%), with the remaining as P-containing heteroatoms (CHONP and CHOSP, 3.0%) and CHOP (3.4%). The much lower occurrence of P-containing compounds was expected with ESI+ mode due to the preferential ionization of N-containing compounds (Ohno et al., 2016; Lu and Liu, 2019). An average of 190 formulas were identified in each sample, with the greatest number in the deep (x̄ = 365 formulas) and shallow groundwater (x̄ = 337 formulas, Table S2 ).

The relative abundance of different compounds in surface water was compared during the pre- and post-landfall and the recovery period of 5 months (i.e., February 2018, Figure 5 , Table S3 ). Pre-landfall surface water was dominated by CHON and CHO compounds. Combined N-containing compounds accounted for ~72% of the compounds detected in surface water prior to landfall. Comparatively, on average, the post-landfall relative abundance of CHO increased by 10%, while that of CHON decreased by ~12%. Also, the average relative abundance of S-containing compounds nearly doubled (10.5%) from pre-landfall (6.6%), with most compounds occurring as NS heteroatoms. After 5 months, the composition of DOM compounds had largely returned to pre-landfall conditions with N-containing compounds once again comprising >70% of the compounds detected, mostly as CHON, and CHO making up most of the rest. However, the relative abundance of S-containing compounds (10.6%) remained like in the post-Harvey samples.

The weighted average molecular element ratios of oxygen to carbon (O/C_w), hydrogen to carbon (H/C_w), and nitrogen to carbon (N/C_w) provide additional insight into DOM composition changes after the landfall. The most oxygenated DOM (i.e., greatest O/C_w) occurred in one February 2018 sample (0.299) and one pre-landfall sample (0.256), while the least oxygenated DOM (i.e., lowest O/C_w, 0.159) occurred on day 10 ( Figure 3D ; Table S2 ). All other surface water samples had relatively consistent O/C_w values (0.214–0.232). These O/C_w values were greater than those previously observed in nearby Nueces Bay during a major flood event (x̄ = 0.185), flood recession (x̄ = 0.176), and more typical low-flow conditions (x̄ = 0.182 and 0.178) (Douglas et al., 2021a). This study had similar surface water unweighted O/C elemental ratios (0.22–0.26) as those observed across four rivers located to the northeast of this study’s location (0.21–0.27), before (at baseflow) and following Hurricane Harvey (at high flow conditions) (Lu and Liu, 2019). The February 2018 peak and day 10 drop in O/C_w coincided with the lowest H/C_w values observed (1.41 and 1.46, respectively). The highest surface water H/C_w values (i.e., the most saturated DOM samples) were observed pre-landfall (>1.60), on days 7, 8, and 11 (1.62–1.69, cluster 2 in post-landfall samples below), and in the other February 2018 sample (1.78). These values are within the upper range of H/C_w values previously observed in the Nueces Bay study at flood peak (1.44–1.77, x̄ = 1.66) and flood recession (1.44–1.84, x̄ = 1.67) (Douglas et al., 2021a) and have similar unweighted elemental H/C ratios (1.48–1.75) to those observed in the monitored rivers at baseflow and high flow conditions (1.46–1.60) (Lu and Liu, 2019). All other post-landfall samples had H/C_w values like the H/C_w values in Nueces Bay during low-flow conditions (x̄ = 1.54 and 1.55). The surface water N/C_w followed similar trends to the H/C_w throughout the study period.

For the groundwater DOM samples, 757 unique compounds were detected in the deeper formation (n = 3), 1,012 compounds in the shallow groundwater (n = 3), and 198 compounds in porewater (n = 2); of these, 611 (81%), 803 (79%), and 153 (77%), respectively, were assigned formulas. Surface water had greater relative abundances of CHON than groundwater samples, while groundwater samples had more CHO, CHOS, and CHOP ( Figure 6 ). Surface water and porewater had similar overall DOM compositions. In contrast, shallow groundwater and deeper groundwater were, on average, composed of a greater proportion of CHO (38%–39%) and S-containing compounds (15%–19%) and had the least N-containing compounds (<50%). Shallow groundwater DOM was generally less oxygenated (O/C_w x̄ = 0.20) and less saturated (H/C_w x̄ = 1.60) than porewater (O/C_w x̄ = 0.23, H/C_w x̄ = 1.67) and DGW DOM (O/C_w x̄ = 0.27, H/C_w x̄ = 1.70, Table S2 ). All groundwater endmembers had generally more saturated DOM than surface water samples, while porewater DOM had similar oxygenation to surface water DOM. Deeper groundwater had the highest percentage of labile compounds (MLB_L x̄ = 82.3%), followed by porewater (MLB_L x̄ = 79.9%) and shallow groundwater (MLB_L x̄ = 70.9%, Table S3 ). Groundwater had generally higher percentages of labile compounds than surface water (MLB_L x̄ = 56.9%).

Cluster analysis and PCA of all surface water and groundwater samples revealed three principal components accounting for ~50% of the variance and comprise 4 clusters that can be explained by differences in groundwater composition (PC1 and PC3) and the DOM composition post-Harvey (PC2, Figures 7A, B ). Principal component 1 accounts for 24.7% of the variance and is driven by the deep groundwater DOM composition (i.e., positive PC1 scores, cluster 1 in Figure 7A ). The scores of PC1 were positively correlated with DOC (p-value < 0.01) and negatively correlated with N/C_w (p-value < 0.001), indicating that deep groundwater had higher concentrations of N-depleted OM. Principal component 2 accounts for 13.4% of the total variance and is driven by the difference in DOM composition immediately following Hurricane Harvey (positive PC2 scores, cluster 2) from pre-landfall and February 2018 surface water, shallow groundwater, and porewater composition (negative PC2 scores, cluster 3, Figures 7A, B ). On the whole, post-landfall surface water samples contained less CHON and more CHO than porewater and pre-landfall and February 2018 surface water samples ( Figures 5 , 6 ). The scores of PC2 positively correlated with pH and S/C_w (p-values < 0.01) and negatively correlated with H/C_w, N/C_w, and NH₄ ⁺ (p-values < 0.05), indicating that post-landfall DOM was less saturated and processed than other samples. The only factor to have a positive relationship with water level (p-value = 0.05) and a weak positive relationship with Oso Creek discharge (p-value = 0.07) was PC2. The difference in shallow groundwater DOM composition was revealed in PC3, which accounts for 12.4% of the variance and is mainly driven by the unique DOM composition of pre-landfall shallow groundwater (cluster 4, positive PC3 scores, Figure 7B ) that had less N (%CHON 30.6, N/C_w 0.06) and a lower percentage of labile compounds (MLB_L 66.3) compared to the other groundwater samples. Principal component 3 was positively correlated with NO₃ ⁻ (p-value < 0.01), weakly correlated with DOC (p-value = 0.06), and negatively correlated with O/C_w and N/C_w (p-values < 0.05).

Cluster analysis and PCA of only post-landfall daily molecular data revealed two principal components and three different groups: cluster 1 (days 4, 5, 6, 9, 10, 12, 22, and 28), cluster 2 (days 7, 8, and 11), and cluster 3 (days 13, 14, and 35) ( Figure 8 ). Explaining ~20% of the variance in the post-landfall DOM composition, PC1 is driven by the difference of days 7, 8, and 11 (positive PC1 scores, cluster 2) from the other post-landfall samples. Cluster 2 had 266 (77%) assigned compounds with a higher presence of CHON and lower CHO than cluster 1 ( Figure 9A ). Samples in PC1/cluster 2 ( Figure 9B ) were positively correlated with H/C_w (p-value < 0.01) and N/C_w (p-value = 0.04) and weakly positively correlated with Oso Creek discharge (p-value = 0.06). In general, cluster 2 had more saturated DOM (H/C_w x̄ = 1.65) with a greater labile nature (MLB_L x̄ = 55.9%, MLB_wL x̄ = 75.4%) than either cluster 1 (H/C_w x̄ = 1.54, MLB_L x̄ = 53.7%, MLB_wL x̄ = 58.8%) or cluster 3 (H/C_w x̄ = 1.56, MLB_L x̄ = 55.2%, MLB_wL x̄ = 61.2%).

An additional ~15% of the DOM variance is explained by PC2 and is driven by the difference between the more immediate post-landfall samples (negative PC2 scores) and the later post-landfall samples (positive PC2 scores), which loosely match cluster 1 and cluster 3, respectively. The PC2 scores were positively correlated with DOC (p-value = 0.03) and water level (p-value < 0.01) and negatively correlated with NO₂ ^– (p-value = 0.05), NH₄ ⁺ (p-value < 0.01), and Oso Creek discharge (p-value = 0.05). Cluster 1 samples had 313 assigned formulas comprised of mainly CHON, CHO, and CHONS ( Figure 9A ).

Cluster 3, with 270 (72%) assigned compounds, had a similar presence of CHON and CHO as cluster 2, but with a higher percentage of CHONS than clusters 1 and 2 ( Figures 6 , 9 ). As mentioned above, samples part of cluster 3 had less saturated DOM with a less labile nature than cluster 2 while slightly more saturated DOM with a more labile nature than cluster 1. The grouping of DOM constituents on the van Krevelen diagram regions shows similarities among the other two clusters, a mix of marine and terrestrial influences ( Figure 9 ).