Massive pollutants released to Galveston Bay during Hurricane Harvey: Understanding their retention and pathway using Lagrangian numerical simulations

a r t i c l e i n f o


Introduction
Pollutant release frequently happens during storm events, especially those accompanied by strong precipitation. The increasing frequency of extreme precipitation events under a warming climate and a more humid atmosphere (Knight and Davis, 2009;Donat et al., 2016;Pfahl et al., 2017) makes the storm-related pollutant release even more threatening to coastal ecosystems in the future. Massive wastewater, nutrient, bacteria, heavy metal, or petrochemical products can be washed away by surface runoff or spilled due to flooding and discharged into receiving waters. Their influence on coastal environments can be catastrophic, particularly for coastal embayment where water exchange with the coastal ocean is slow and pollutants can stay for a long time. Damages to the water quality, marine environment, marine mammals, and fishery due to released pollutants have been extensively observed (Weyhenmeyer et al., 2004;Cardoso et al., 2008;Wetz and Yoskowitz, 2013). Recovery time from such extreme events for the receiving waters in terms of hydrodynamics or ecosystem health can take months or even years, depending on the amount of freshwater load, pollutant concentration, flushing capacity, and resiliency of the ecosystem (Paerl et al., 2001).
Environmental assessment for pollutant release usually focuses more on the total pollutant load while rarely taking into account the timing of release. Often, more pollutants are released during storm discharge, but the flushing is also stronger during the storm. Their influence on the coastal water quality is therefore not necessarily linearly proportional to the total loading (Taylor et al., 2011). As suggested by Dettmann (2001), more fraction of pollutants will be exported out of coastal systems when the flushing capacity is higher. Taking the massive pollutants released from the Houston metropolitan area in Texas to the adjacent Galveston Bay during Hurricane Harvey as an example, we show here that the release timing is critically important and the susceptibility of coastal waters to pollutant can be more serious than expected when pollutants are released after the storm discharge.
Hurricane Harvey, the wettest tropical cyclone on record in the U.S., made landfall on August 26, 2017 along the mid-Texas coast as a Category 4 hurricane and brought unprecedented rainfall to the Texas-Louisiana coast, with a return period of the peak 3-day precipitation exceeding 1000 years (van Oldenborgh et al., 2018). Intense rain with the daily precipitation averaged over the bay area larger than 50 mm lasted for 5 days (i.e., August 26-30). Maximum accumulative precipitation reached 1539 mm (60.58 in.) (Mathews, 2019), causing more than 80 deaths and over 150 billion dollars of economic loss. It was estimated that Harvey delivered 14 Â 10 9 m 3 of freshwater (~3.7 times of the bay's volume) and deposited 9.9 Â 10 7 metric tons of sediment (equivalent to 18 years of average annual sediment load) to Galveston Bay (Du et al., 2019a,b). The bay became virtually fresh for a few days, and salinity recovery inside the bay took about 2 months on average .
To make things worse, many petrochemical facilities were flooded, resulting in chemical pollutant leak or release (Fig. 1a). The flooding and the subsequent pollutant release are of great concern, since Houston is known as the second-largest petrochemical industry hubs in the world (Santschi et al., 2001) and the fourth largest city in the U.S. in terms of population size. Harvey was estimated to cause release of 0.57 Â 10 6 tons of raw sewage (Phillips, 2018) and more than 22,000 barrels of oil, refined fuels and chemicals (Flitter and Valdmanis, 2017) to Galveston Bay. Harvey's aftermath lasted for a long time, e.g., drastic mortality and slow recovery of oysters (Christine Jensen, personal communication) and excessive skin problems for dolphins (Stuckey, 2017).
To date, some questions regarding the impacts of released pollutants during Harvey are still not answered. For example, how were the released pollutants dispersed inside and outside the bay? How long did they stay inside the bay and where did the pollutants aggregate? What are their pathway differences between the normal condition and during storm discharge? Understanding these questions is essential for environmental assessment, water quality management, and ecosystem restoration.
In this study, we used a Lagrangian particle-tracking method coupled with a validated 3D hydrodynamic model to examine the retention and pathway of pollutants released during Hurricane Harvey. Due to the random nature of particle movement and the large number of released particles, it is unpractical to analyze the pathway for each particle. To this end, we introduce a new transport timescale, called local exposure time, to describe the spatially varying susceptibility and to synthesize the mean characteristics of particle dispersion. Considering the increasing intensity and frequency of precipitation events, this study for Hurricane Harvey will be instructive for future research.
The paper is structured as follows: Section 2 briefly introduces the 3D numerical model, the coupled Lagrangian particle-tracking method, and two transport timescales used to quantify the particle retention and describe their pathway. Section 3 presents the results of Lagrangian simulations, with special focus on the difference between post-storm release and storm release. Section 4 discusses the importance of timing on the particles' retention and pathway, as well as the role of shelf and ocean circulations for the dispersion of particles after they exited the coastal system.

Hydrodynamic model
We employed the Semi-implicit Cross-scale Hydroscience Integrated System Model (SCHISM: Zhang et al., 2015Zhang et al., , 2016, an opensource community-supported modeling system based on unstructured grids, derived from the early SELFE model (Zhang and Baptista, 2008). SCHISM uses a highly efficient semi-implicit finite-element/finite-volume method with a Eulerian-Lagrangian algorithm to solve the turbulence-averaged Navier-Stokes equations, under the hydrostatic approximation. It uses the generic length-scale model of Umlauf and Burchard (2003) with the stability function of Kantha and Clayson (1994) for turbulence closure. One of the major advantages of the model is that it has the capability of employing a very flexible vertical grid system, robustly and faithfully resolving the complex topography in estuarine and oceanic systems without any smoothing (Zhang et al., 2016;Stanev et al., 2017;Du et al., 2018a;Ye et al., 2018).
The model domain (Fig. 1c) covers the entire Texas, Louisiana, Mississippi, and Alabama coasts, including the shelf as well as major estuaries (e.g., Galveston Bay). The model grid has a resolution ranging from 40 m in the narrow ship channel of Galveston Bay to 2.5 km on average over the shelf and 10 km in the open ocean. Vertically, a hybrid s-z grid is used, with 10 sigma layers for depths less than 20 m and another 30 z-layers for depths from 20 to 4000 m with shaved cells near the bottom. The bathymetry used in the model is based on the coastal relief model (3 arc-second resolution: https://www. ngdc.noaa.gov/mgg/coastal/crm.html, last access: September 25, 2019). The local bathymetry in Galveston Bay is augmented with 10-m resolution DEM (digital elevation model) bathymetric data (https://catalog.data.gov/dataset/galveston-texas-coastal-digitalelevation-model, last access: September 25, 2019) to resolve the narrow ship channel (150 m wide, 10-15 m deep) that extends from the bay entrance all the way to Port of Houston. When forced by realistic boundary conditions, including the open boundary conditions from FES2014 global tide (Carrere et al., 2015) and global HYCOM model output (https://www.hycom.org/data/glbu0pt08, last access: September 25, 2019), atmospheric forcing from the European Centre for Medium-Range Weather Forecasts (ECMWF: https://www. ecmwf.int, last access: September 25, 2019), and freshwater discharges for 15 rivers, the model gives a good reproduction of the observed hydrodynamic conditions in 2007-2008 inside the Galveston Bay and over the Texas-Louisiana shelf in terms of water level, salinity, temperature, stratification, and shelf currents (Du et al., 2019c).
The model has been applied to simulate the hydrodynamic conditions during Hurricane Harvey, and it reproduced well the dramatic estuarine responses, including the long-lasting elevated water level, extraordinarily strong along-channel velocity, sharp decreases and long recovery of salinity, and huge river plumes on the shelf . The validated hydrodynamic model provides reliable hydrodynamic fields, with which the following Lagrangian simulations are coupled.

Lagrangian particle tracking
A Lagrangian particle tracking method coupled with the 3D hydrodynamic model outputs was used to simulate the dispersion of pollutants. At the junction between San Jacinto River and Buffalo Bayou (Fig. 1d), 1378 particles (neutrally buoyant) were released every day at 00:00 from August 1st to October 1st, 2017. The release location was selected because the most serious pollutant release was from the petrochemical facilities along the Buffalo Bayou and San Jacinto River (Fig. 1a). A random walk was implemented in the particle tracking module to include the influence of the diffusion processes. The movements of particles are governed by advective and diffusive transport processes.
where (X, Y, Z) is the location of the particle; U, V, and W are the water velocity components in the Cartesian coordinates of x, y, and z, respectively; n and n+1 indicate the current and next time steps, respectively; Dt is the time interval; R is a uniform random number between À 1 and 1; and K x , K y , and K z are turbulent diffusion coefficients in the x, y, and z directions, respectively.

Transport timescales
Two transport timescales were calculated to quantify the retention of the particles. One is the transit time, which measures the duration of a particle staying inside a defined domain and is calculated as the time difference between entering and exiting the domain (Shen and Haas, 2004). The mean transit time (/Þ averaged over the N particles is calculated as where t 2 is the time when the particle exits the domain (i.e., Galveston Bay in this study) for the first time, and t 1 is the time when the particle enters the domain. The transit time here does not consider the returning of particles. The local exposure time (LET) is the other timescale used. It is a new concept derived from the traditional exposure time that measures the overall lifetime a particle spends inside a given domain. The exposure time includes the duration after the particle returns into the domain (Delhez, 2006). LET is the mean exposure time of a set of particles within a defined region. Different from traditional exposure time that gives one scalar value, LET allows us to examine the spatial variability as we can separate the domain into many small regions. LET can be calculated as, where r i = 1 when the particle i is inside the defined region and r i = 0 when the particle i is outside the defined region. The integration window with respect to t is from the particle release time to the end of model run. For this study, the hydrodynamic output from the numerical model covers the period of July 1 to December 31, 2017 . Another practical way to calculate the LET from numerical results is to register the time when a particle exits or enters the defined region, integrate the duration of all the visits, and average for all particles (Fig. 2). In this study, we divided the entire domain (including the bay and adjacent shelf) into 1 km Â 1 km square regions and calculated the LET for each region. Both particle's pathway and retention affect the degree to which released pollutants would influence local water quality and ecosystem health. Since it is impractical to examine the pathway of all released particles, LET provides a succinct and quantitative description of all particles. LET is the combined result of the flushing/exchange efficiency and the possibility of particles reaching the given region. It is, therefore, suitable to use LET as a measure of the susceptibleness of any given region to the released pollutant. A longer LET indicates the region is more susceptible to the pollutant. For a given region, a longer LET can be caused by several reasons including: (a) more particles passing through; (b) more visits by each particle (e.g., back-forth moving due to tide); (c) longer retention time for each visit (e.g., due to slow current).

Particles released during storm
Particles released at the beginning of Harvey discharge (i.e., August 27, 2017) were quickly flushed out of the bay, with a median transit time of 1.5 days and 90% of the particles exiting the bay within 2 days (Fig. 3a-h). Particles moved seaward along the longitudinal axis of the bay without tidal (back-and-forth) movement (Fig. 4c). After exiting the bay, due to strong seaward momentum, particles moved offshore, reaching as far as 50 km off the bay entrance, consistent with remotely sensed sediment plume . The fast seaward movement was primarily caused by the strong ebbing along-channel current during Harvey. The velocity measured at the bay entrance and the upper bay showed the current (tidal + subtidal) was in seaward direction during the entire storm discharge period, with the maximum ebbing velocity exceeding 3 m s À1 in the upper bay (Du et al., 2019a).
As a result, the transit time was rather small during the storm discharge period, with a minimum value of 1 day. The transit time shows a linearly decreasing trend from August 1 to August 30 (Fig. 4a), because particles released before Harvey were subject to quick flushing during the storm discharge. Therefore, the transit time for particles released before Harvey is merely the time difference between the particle release time and the beginning of storm discharge. It suggests that an episodic flooding event can efficiently refresh the entire bay, thus playing an important role in the overall water renewal.

Particles released after storm
The movement, pathway, and transit time of particles released after the storm discharge were distinctly different from those released at the beginning of Harvey. Particles moved seaward slowly, with a median transit time of 60-90 days (Fig. 3i-p). Particles moved back and forth under the influence of tidal currents and were dispersed over the entire bay (Fig. 4e). It took a long time for particles to exit low-flushing regions such as Trinity Bay and East Bay. After exiting the bay, particles tended to move downcoast and mostly concentrated on the inner shelf (Figs. 2k-o and 3d). Surprisingly, the transit time quickly reverted to that under normal condition (~80 days) right after the storm discharge was termi- nated. It suggests a quick recovery in flushing or water exchange after the storm, despite the fact that it took about two months for the salinity inside the bay to recover to the pre-storm condition .

LET to quantify the susceptibleness
For the particles released during storm (i.e., August 26-30), LET was small (Fig. 5a). With a maximum value of about 0.3 h, LET was less than 0.1 h for most regions, suggesting a fast movement and short retention of particles. LET was larger along the ship channel and decreased as moving away from the channel, indicating that most particles moved seaward along the ship channel. It is necessary to point out that stronger along-channel velocity made more particles move along the ship channel, resulting in more exposure time of particles and thus larger LETs near the ship channel.
For the particles released during September 3-7 (i.e., after the storm), LETs were at least one-order larger than those of storm release (Fig. 5b). LET reached a maximum value of about 6 h in the upper bay and had large values in Trinity Bay. It suggests that particles tended to move into Trinity Bay and stayed there for a long time, which was related to the small tidal range and slow water renewal in Trinity Bay. Salinity measurement and numerical modeling have confirmed the negligible tidal signal and slow recovery in salinity inside the Trinity Bay (Du et al., 2019a;. LET became smaller toward downstream, primarily due to a strong tidal exchange between the lower bay and the shelf. Particles exiting the bay would have less chance to return due to shelf transport, very different from that in the middleupper bay, where particles moved back and forth with tidal cycles.

Importance of pollutant release timing
The model results draw our attention to an important but previously not well-recognized aspect concerning the pollutant susceptibility, that is, the release timing matters. Environmental assessments typically focused on the total pollutant loading (Brezonik and Stadelmann, 2002;Nazahiyah et al., 2007;Tiefenthaler et al., 2008) and paid little attention to the release timing. We demonstrate here that pollutants released after the storm will be more influential on the water quality and ecosystem health, as they will stay much longer inside the bay, than those released during storm (Fig. 5) although the amount of pollutants released after the storm is usually much smaller than that during storm. It should be noted that the loading after Harvey was still significant due to a large population and dense petrochemical industries around Galveston Bay. Furthermore, several reservoirs were controlled to release polluted water slowly. For example, Barker Reservoir, located west of Houston, released freshwater for over 40 days after the storm (Du et al., 2019a). We can imagine that pollutants in the reservoir waters were also released slowly and steadily, which, combined with those from the city and groundwater along the coastline, might deteriorate the bay's water quality well after the storm discharge.
The underlying mechanism responsible for the timing sensitivity lies in the overall ocean-estuary exchange or the flushing capacity. Slow flushing of the bay is the primary factor amplifying the timing issue. The storm discharge and the resulting seaward outflow continued for only 5-7 days and then the back-and-forth tidal current resumed. Pollutants released after the storm discharge is  subject to long retention inside the bay. If the flushing time of a system under normal condition is small, the timing sensitivity can be of less importance. For many estuarine systems in the Gulf of Mexico, however, flushing times are relatively long due to small tidal range and narrow outlets (Rayson et al., 2016;Du et al., 2018b). Relatively long flushing times are also expected for small coastal embayments with small tidal range and little freshwater discharge.
It is worthy to note estuarine systems with relatively wide mouth(s), strong tide, and large volume are unlikely to have such dramatic influence. For example, in Chesapeake Bay or San Francisco Bay, two of the largest estuaries in the U.S., the flushing time varies in a less dramatic manner, primarily due to the large volumes, and it is unlikely that the entire bay will be flushed out during a storm event (Walters et al., 1985;Du and Shen, 2016). As a result, the impact of the pollutants released before or after the storm discharge will be similar, and the amount of pollutants (loading) will become more important, although the release timing can affect some local regions, particularly near the release locations.
The fate and pathway of released pollutants in coastal bays are quite different from that in river systems where released pollutants are generally controlled by diffusion and one-way advection (Whitehead et al., 1986;Chapra and Whitehead, 2009). Complex geometric and bathymetric features in coastal bays, together with the barotropic and baroclinic interaction with shelf oceans, make the pathway and susceptibility of pollutants difficult to predict and quantify if without numerical tools and useful indexes. Using the timescale LET would significantly simplify the environmental assessment. One application of the LET is to examine the different fate and pathways of pollutants released at different locations. Numerical experiments show that LET varies greatly depending on the release location (Fig. 6). For instance, pollutants released at the Trinity River mouth tended to have more influence on East Bay (Fig. 6a), while pollutants released at Clear Lake (i.e., western shore) tended to have more impact on West Bay (Fig. 6b). It is necessary to note that LET is normalized by the number of particles released and thus indicates the average exposure for unit amount of pollutants. When assessing the susceptibility of a defined region, information about the amount of pollutants should be obtained first.

Importance of shelf and ocean circulation for pollutant dispersion
After exiting the bay, dispersion of the pollutants was directly regulated by the shelf circulation, ocean currents, and the interaction between the shallow shelf and deep ocean. Meso-scale eddies in the ocean can be influential in altering the water exchange between ocean and shelf.
Physical oceanography in the Gulf of Mexico is characterized with intrusions of Loop Currents (Oey et al., 2005), mesoscale eddies derived from Loop Currents (Barkan et al., 2017), and one of the largest river systems in the world, Mississippi River (Rabalais et al., 2002). Despite the expectation that a strong hurricane like Harvey might affect the shelf and ocean circulation, analysis of the satellite-data based sea surface height (SSH) and geostrophic currents shows little change was induced by Harvey (Fig. 7). The SSH saw a marked increase in the coastal waters near Galveston Bay from August 24 to 31, 2017 (Fig. 7d), which was believed to be caused by the addition of large freshwater during Harvey. This analysis suggests, despite the great influence of Harvey on local coastal systems, its influence on the overall ocean circulation was negligible. However, it might be different for other hurricanes. From drifter data and numerical modeling, Curcic et al. (2016) showed Hurricane Issac in 2012 caused significant stokes drift and ocean waves. Oey et al. (2006) showed that a significant warming of Loop Current was induced by Hurricane Wilma in 2005.
After exiting the local estuarine system, pollutant dispersion in the coastal ocean will follow the usually disturbed coastal circulation and their fate may differ from that under normal conditions. It is necessary to point out that the shelf circulation on the Texas-Louisiana shelf has clear seasonality, with downcoast shelf current most of the time except during summer (Cochrane and Kelly, 1986;Cho et al., 1998). Particles released under normal conditions (taking 2007-2008 as an example) were mostly transported downcoast except during July and August (Fig. 8). One interesting pattern in the particle distribution is that the eddies along the shelf break play a role in augmenting the water exchange between the shelf and deep ocean. The mesoscale eddies (length scale~100 km) may persist months in the Gulf of Mexico. For instance, the warm-core rings detach from the Loop Current episodically at an interval of 4-17 months, move slowly westward, and have lifetimes from months to a year (Sturges and Leben, 2000).
In summary, pollutant dispersion is subject to the influence of local bathymetric and geometric features and to the regulation of shelf-ocean circulations. One thing we did not include in this study is the settling or buoyancy property of pollutants. Certain types of pollutants tend to float at the sea surface due to smaller density and hydrophobic nature, while others tend to settle because of heavier density or attachment to suspended sediment. To include these kinetic processes of pollutants, one has to prove the numerical model is reasonably accurate in simulating not only the hydrodynamics but also sedimentary processes. Such efforts should be conducted in future research to facilitate the assessment of certain pollutants. As for the concern of this study, the analysis of passive particles is sufficient to support the major conclusions. Even though the retention time might be moderately or dramatically changed if settling or floating is considered, the pathway and distinct difference due to release timing will very likely maintain the same pattern.

Conclusions
Given the increasing frequency of extreme precipitation events as projected in several recent studies (e.g., Knight and Davis, 2009;Donat et al., 2016;Pfahl et al., 2017), intense pollutant releases are expected to occur more frequently in the future. This study uses a Lagrangian particle-tracking method to examine how massive pollutants released during Hurricane Harvey were dispersed inside the Galveston Bay and over the adjacent shelf. We found distinctively different retention and pathway between particles released during and after the storm discharge. Using LET, it is found that the susceptibility to the released pollutants was at least oneorder greater for post-storm release than storm release. It suggests the timing of pollutant release can be more critical than the total amount of pollutants, particularly for slow flushing, smallvolume estuarine systems. Our study suggests that the environmental agencies pay attention not only to the total pollutant load but also to the timing of pollutant release.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.