Water Circulation Driven by Cold Fronts in the Wax Lake Delta (Louisiana, USA)

Abstract

Atmospheric cold fronts can periodically generate storm surges and affect sediment transport in the Northern Gulf of Mexico (NGOM). In this paper, we evaluate water circulation spatiotemporal patterns induced by six atmospheric cold front events in the Wax Lake Delta (WLD) in coastal Louisiana using the 3-D hydrodynamic model ECOM-si. Model simulations show that channelized and inter-distributary water flow is significantly impacted by cold fronts. Water volume transport throughout the deltaic channel network is not just constrained to the main channels but also occurs laterally across channels accounting for about a quarter of the total flow. Results show that a significant landward flow occurs across the delta prior to the frontal passage, resulting in a positive storm surge on the coast. The along-channel current velocity dominates while cross-channel water transport occurs at the southwest lobe during the post-frontal stage. Depending on local weather conditions, the cold-front-induced flushing event lasts for 1.7 to 7 days and can flush 32–76% of the total water mass out of the system, a greater range of variability than previous reports. The magnitude of water flushed out of the system is not necessarily dependent on the duration of the frontal events. An energy partitioning analysis shows that the relative importance of subtidal energy (10–45% of the total) and tidal energy (20–70%) varies substantially from station to station and is linked to the weather impact. It is important to note that within the WLD region, the weather-induced subtidal energy (46–66% of the total) is much greater than the diurnal tidal energy (13–25% of the total). The wind associated with cold fronts in winter is the main factor controlling water circulation in the WLD and is a major driver in the spatial configuration of the channel network and delta progradation rates.

1. Introduction

Atmospheric cold fronts have an average interval of 3–7 days in the late fall to spring and dominate regional weather along coastal Louisiana [1]. According to a statistical analysis of 40 years of weather data, there are an average of 41 ± 5 cold fronts passing through the area each year between mid-October and the following April [2]. These atmospheric events generally move from northwest to southeast, mostly transiting temperate and subtropical regions, and are characterized by a sequence of spatiotemporal changes in wind speed and direction, barometric pressure, air temperature, and humidity. The associated wind drives the water circulation in estuarine and shelf waters [3,4]. For instance, winds from the southern quadrants prior to a cold front passage tend to drive saltwater intrusion and associated water circulation, particularly through narrow inlets across different estuarine systems in coastal Louisiana (e.g., Calcasieu Lake [5], Lake Pontchartrain [6,7,8], Port Fourchon [9], and Barataria Bay [10]). This wind pattern can also cause either erosion or deposition (e.g., Chenier Plain [11]) depending on (1) the orientation, (2) direction of propagation, (3) angle of approach (the angle to that of the coastline orientation), (4) propagation speed (controlling the duration of wind forcing), and (5) wind field, of the cold front.

One of the areas significantly impacted by cold fronts is the Wax Lake Delta (WLD), which is part of the Atchafalaya-Wax Lake complex located in the eastern region of the largest delta plain in the Northern Gulf of Mexico (NGOM). This coastal region is strongly influenced by cold fronts as well as tropical cyclones (e.g., [3,12,13]). The WLD is a river-dominated delta [14,15,16] initially formed by sediment deposition after a flood-control canal was dredged in 1942 (i.e., Wax Lake Outlet [17]) to regulate excess water flow in the Mississippi River during the spring peak discharge, when it becomes a flooding threat to coastal urban areas. Overall, the outlet diverts approximately 30% of the flow and sediment from the Atchafalaya River directly down to the Atchafalaya Bay (Figure 1).

The WLD has been expanding as diverted sediment fills the river mouth area [18] and successionally builds hydrologically connected delta lobes throughout a network of channels and bars [14,19,20,21]. Owing to the differences in sediment deposition rate and delta lobe age, these lobes are colonized by wetland vegetation in different stages of succession and species composition [22,23] that, in turn, influence sediment deposition dynamics [21,24,25]. Given the variable influence of the Mississippi River flow, the Atchafalaya-WLD complex sediment composition is dominated by sand (~70%). Wave actions, water level fluctuations, and winds associated with atmospheric cold fronts are responsible for the resuspension of fine-grained sediments and their subsequent transport to the adjacent continental shelf. Both the Atchafalaya and Wax Lake deltas continue to prograde by filling extensive areas across the south-central region of the Atchafalaya Bay as they advance into the continental shelf [26,27]. The rapid evolution of the growing delta complex represents one of Louisiana’s most dynamic coastal environments and serves as an example of the potential and desired outcomes in river diversion projects to promote land-building in the lower Mississippi River [28].

Further, the current WLD, with vegetation development, is an outstanding example of how thriving wetland communities can recover from the impacts of tropical cyclones and strong seasonal cold fronts. These severe weather events can produce excess flooding during high storm surges and subsequent drainage of the system that controls delta morphological evolution and wetland structural and functional properties [11,29]. Within this deltaic region, distributary channels bounded by the vegetated lobes or islands, may be submerged or emergent depending on the Mississippi River stage, wind (magnitude/direction), and the local wetland hydroperiod regime (frequency/duration of inundation and water depth [24]). During high flooding periods caused by seasonal variations or weather events, the natural levees can be overtopped thus increasing water exchange between the channels and wetland habitats inside the islands [27,30,31,32,33]. Hence, measuring sediment and nutrient transport across these inter-distributary channels requires an assessment of changes not only in bathymetry, water flow, and spatiotemporal shifts in vegetation density/thickness (i.e., flow resistance) [24,33,34,35,36,37], but also an understanding of the seasonal pulsing impact of winter cold fronts regulating flow distribution and circulation inside the WLD.

Because winter cold fronts have a significant impact on the resuspension and transport of materials in the WLD region [38], it is critical to understand water circulation during these events to determine the net sediment, nutrient and carbon transports within the WLD proper [36,37,39,40] and between the Atchafalaya-WLD system and the adjacent shelf area [38,40,41,42,43]. However, direct hydrological measurements are particularly challenging during severe weather such as cold fronts. Thus, high-resolution numerical modeling of hydrodynamic processes and the use of (albeit limited) field data obtained before and after cold fronts can help discern the effects and mechanisms controlling flooding frequency and duration, water level variability, and the transport of materials across boundaries.

In view of the above discussion and the need for a better understanding, we perform a set of numerical model experiments to investigate the key physical mechanisms responsible for water circulation and transport in the Atchafalaya-WLD system. We use a modified version of the semi-implicit Estuary and Coastal Ocean Model (ECOM-si [44]). The numerical experiments are designed to (1) evaluate the relative importance of distributary channels in diverting and attenuating water flow, (2) determine the differences in water level before and after cold frontal passages, (3) evaluate how the interaction between the tidal regime and the water level induced by cold fronts controls water flushing rates from the WLD, and (4) partition the relative energy contribution among different water circulation drivers. In this work, we ignore the overall river discharge effect on water level change inside the delta since the first order variability of river discharge is seasonal in terms of its magnitude, i.e., the impact occurs over a much longer time scale than that of individual cold front events; therefore, most of the cold front induced circulation and water level fluctuations are wind-dominated.

2. Data and Methodology

This study uses a numerical experiment to evaluate how the WLD hydrodynamics varies under cold front weather conditions. Both model calibration and validation were performed followed by model simulations that include several cold fronts and the analysis of water flow fields. Energy partitioning was then applied to determine the relative importance of different water flow components.

2.1. Data

The data used in model calibration and parametrization include local and regional measurements of meteorology, hydrodynamics, bathymetry, horizontal/vertical velocity profiles, river discharge, and LIDAR topography.

Meteorological and current velocity data were obtained from three offshore stations of the Wave-Current-Surge Information System (WAVCIS, http://www.wavcis.lsu.edu/, accessed on 7 March 2022) located at water depths of 5 (CSI3), 20 (CSI6), and 15 (CSI9) m along the Louisiana shelf (Figure 1). Additional data from two deployed instruments measuring current velocity and water level were also used; the first deployment (Delta #1) was placed in the fourth channel of WLD located east of the delta (water depth is 1.7 m) while the second one was deployed in the Big Hog Bayou (water depth, ~1.4 m deep, Figure 1) between the WLD and the Atchafalaya Delta (AD). Most of the WLD region has very shallow water on the order of 1–2 m, except in the channels. The central channel has a maximum depth of ~18 m due to erosion from strong river flows. The water depth shallows quickly toward the rim of the delta. The fan-shaped WLD has a length of about 11 km along its axis and 13 km in the lateral direction. The channel width ranges from 0.25 km to 2 km (Figure 2).

The WLD bathymetry data for model implementation were obtained from two sources: (1) the National Geophysical Data Center’s (NGDC) 3 arc-second (~90 m) U.S. Coastal Relief Model (CRM: http://www.ngdc.noaa.gov/mgg/coastal/crm.html, accessed on 7 March 2022); which integrates offshore bathymetry with land topography to provide a comprehensive and seamless representation of the U.S. coastal zone, and (2) high-resolution bathymetry surveys obtained using a small research vessel in 2011, 2012, 2013 (between longitudes 91.5848° W and 91.292° W; and latitudes 29.3647° N and 29.6466° N, [45]); these in situ high-resolution bathymetric measurements helped to improve data accuracy, especially in the channels.

The daily mean river discharge data were obtained from the USGS stations # 07381590 (Wax Lake Outlet at Calumet, LA) and #07381600 (the Lower Atchafalaya River at Morgan City, LA. Figure 1); the latter is located near the Atchafalaya Bay River mouth (

Table 1. Meteorological and oceanographic data from different stations (Figure 1). (* indicates available data).

Station	Location		Data
	Latitude	Longitude	Wind	Barometric Pressure	Air Temperature	Water Level	Current
Delta #1	29.506357°	−91.472336°			*	*	*
Big Hogs Bayou	29.518045°	−91.354841°			*	*	*
Small Tripod	29.501126°	−91.477912°				*
Large Tripod	29.517809°	−91.354793°				*
CSI3	29.441°	−92.061°	*	*	*	*	*
CSI6	28.867°	−90.483°	*	*	*	*	*
CSI9	29.1015°	−89.9782°	*	*	*	*
LAWMA	29.448333°	−91.336667°	*	*	*	*
BE	29.666667°	−91.236667°	*	*	*	*
FRWL	29.555°	−92.305°	*	*	*	*
GI	29.263333°	−89.956667°	*	*	*	*
PS	29.178333°	−89.258333°	*	*	*	*
PI	28.931667°	−89.406667°	*	*	*	*
Port Fourchon	29.113333°	−90.198333°				*
WL	29.697778°	−91.372778°				*
LAR	29.692611°	−91.211833°				*

).

2.2. Methodology

The numerical model experiments use a fully nonlinear, three-dimensional, curvilinear coordinate, finite-difference, and structured-grid numerical model based on shallow water (hydrostatic) hydrodynamics. This model is a modified version [46,47] of the ECOM-si initially developed based on a version of the Princeton Ocean Model (POM [48]). The earlier version of this model used an orthogonal curvilinear coordinate system (e.g., [49]) while the modified version by [46], extended it to a non-orthogonal curvilinear coordinate system, allowing more flexibility than the earlier version. The advantage and limitation of the model were discussed in [47] and [50]. Several previous applications have demonstrated the model’s capability in reproducing estuarine water circulation, freshwater discharge plumes, water exchange, suspended sediment transport, and water quality [49,50,51,52]. The governing equations [46] for momentum, continuity, water temperature, and density are defined in a horizontally non-orthogonal, and vertically stretched sigma-coordinate system:

$$\begin{gathered} \frac{\partial DJ{u}_1}{\partial t}+\frac{\partial DJ\widehat{U}{u}_1}{\partial \xi }+\frac{\partial DJ\widehat{V}{u}_1}{\partial \eta }+\frac{\partial J\omega u_1}{\partial \sigma }-D{h}_2\widehat{V}\left[v_1\frac{\partial }{\partial \xi }\left(\frac{J}{h_1}\right)-u_1\frac{\partial }{\partial \eta }\left(\frac{J}{h_2}\right)+Jf\right] \\ -D{h}_2{u}_1{v}_1\frac{\partial }{\partial \xi }\left(\frac{h_3}{h_1{h}_2}\right) \\ =-h_{2}gD\frac{\partial \zeta }{\partial \xi }+\frac{g{h}_{2}D}{{\rho }_o}\frac{\partial D}{\partial \xi }{{\displaystyle \int }}_{\sigma }^0\sigma \frac{\partial \rho }{\partial \sigma }d\sigma -\frac{g{h}_2{D}^2}{{\rho }_o}\frac{\partial }{\partial \xi }{{\displaystyle \int }}_{\sigma }^0\rho d\sigma +\frac{1}{D}\frac{\partial }{\partial \sigma }\left(K_m\frac{\partial J{u}_1}{\partial \sigma }\right)+DJ{F}_x \end{gathered}$$

(1)

$$\begin{gathered} \frac{\partial DJ{v}_1}{\partial t}+\frac{\partial DJ\widehat{U}{v}_1}{\partial \xi }+\frac{\partial DJ\widehat{V}{v}_1}{\partial \eta }+\frac{\partial J\omega v_1}{\partial \sigma }+D{h}_1\widehat{U}\left[v_1\frac{\partial }{\partial \xi }\left(\frac{J}{h_1}\right)-u_1\frac{\partial }{\partial \eta }\left(\frac{J}{h_2}\right)+Jf\right] \\ -D{h}_1{u}_1{v}_1\frac{\partial }{\partial \eta }\left(\frac{h_3}{h_1{h}_2}\right) \\ =-h_{1}gD\frac{\partial \zeta }{\partial \xi }+\frac{g{h}_{1}D}{{\rho }_o}\frac{\partial D}{\partial \eta }{{\displaystyle \int }}_{\sigma }^0\sigma \frac{\partial \rho }{\partial \sigma }d\sigma -\frac{g{h}_1{D}^2}{{\rho }_o}\frac{\partial }{\partial \eta }{{\displaystyle \int }}_{\sigma }^0\rho d\sigma +\frac{1}{D}\frac{\partial }{\partial \sigma }\left(K_m\frac{\partial J{v}_1}{\partial \sigma }\right)+DJ{F}_y \end{gathered}$$

(2)

$$\frac{\partial \zeta }{\partial t}+\frac{1}{J}\left[\frac{\partial }{\partial \xi }\left(DJ\widehat{U}\right)+\frac{\partial }{\partial \eta }\left(DJ\widehat{V}\right)\right]+\frac{\partial \omega }{\partial \sigma }=0$$

(3)

where

$$\omega =w-\sigma \left(\widehat{U}\frac{\partial D}{\partial \xi }+\widehat{V}\frac{\partial D}{\partial \eta }\right)-\left[\left(1+\sigma \right)\frac{\partial \zeta }{\partial t}+\widehat{U}\frac{\partial \zeta }{\partial \xi }+\widehat{V}\frac{\partial \zeta }{\partial \eta }\right],$$

(4)

and $\xi , \eta , and \sigma$ are defined as

$$\begin{array}{ccc}\xi =\xi \left(x, y\right),& \eta =\eta \left(x, y\right),& \sigma =\frac{z-\zeta }{H+\zeta }.\end{array}$$

(5)

Here $\sigma$ varies from −1 at $z=-H$ to 0 at $z=\zeta$; $x, y, \mathrm{and} z$ are the eastward, northward, and upward axes of the orthogonal Cartesian coordinates; $\zeta$ is the surface elevation; $H$ is the still water depth; $D$ is the total water depth ($H+\zeta$); f is the Coriolis parameter; $g$ is the gravitational acceleration, and $K_m$ is the vertical eddy viscosity coefficient. The coefficient $K_m$ is calculated using the modified Mellor and Yamada level 2.5 turbulent closure scheme [53,54]. $F_u \mathrm{and} F_v$ represent the horizontal momentum diffusion terms calculated by Smagorinsky’s formula [55] where the horizontal diffusion is proportional to the product of horizontal grid sizes.

In this non-orthogonal curvilinear coordinate system, $u_1 \mathrm{and} v_1$ are the $\xi \mathrm{and} \eta$ components of the velocity that can be converted back to the $x \mathrm{and} y$ components ($u \mathrm{and} v$) of the velocity using the following [46]:

$$\begin{array}{cc}{u}_1=\frac{h_2}{J}\left(x_{\xi }u+y_{\xi }v\right),& v_1=\frac{h_1}{J}\left(x_{\eta }u+y_{\eta }v\right),\end{array}$$

(6)

where $J$ is the Jacobian function $J=x_{\xi }{y}_{\eta }-x_{\eta }{y}_{\xi }$. The subscripts $\xi \mathrm{and} \eta$ indicate partial derivatives with respect to $\xi \mathrm{and} \eta$, respectively. The factors $h_1 \mathrm{and} h_2$ of the coordinate transformation are defined as

$$\begin{array}{cc}{h}_1=\sqrt{x_{\xi }^2+y_{\xi }^2},& h_2=\sqrt{x_{\eta }^2+y_{\eta }^2},\end{array}$$

(7)

and $\widehat{U} \mathrm{and} \widehat{V}$ are given as

$$\begin{array}{c}\widehat{U}=\frac{1}{J}\left(h_2{u}_1-\frac{h_3}{h_1}{v}_1\right),\\ \widehat{V}=\frac{1}{J}\left(h_1{v}_1-\frac{h_3}{h_2}{u}_1\right), \end{array}$$

(8)

where $h_3=y_{\xi }{y}_{\eta }+x_{\xi }{x}_{\eta }$. $\rho \mathrm{and} {\rho }_o$ are the perturbation and reference densities, which satisfy ${\rho }_{total}=\rho +{\rho }_o$. The momentum equations are nonlinear with a variable Coriolis parameter. Since the storm surge generated by strong winds is mainly barotropic, the temperature and salinity equations are not included in this study.

3. Model Implementation

3.1. Setup

The model domain encompasses the Louisiana inner shelf from the Texas-Louisiana border to the Louisiana Bight and Mississippi River bird-foot delta (Figure 3). It covers the Atchafalaya Bay system which encompasses five contiguous bays: Vermillion Bay, West and East Cote Blanche Bays, Atchafalaya Bay, and Fourleague Bay [56,57]. The model domain is fan-shaped and centered at the WLD including the coastal, shelf, and shelf break regions. The model is divided into $380\times 250$ cells in the curvilinear horizontal orthogonal coordinate system (Figure 3$;-{93.68}^{°} \mathrm{W},-{89.21}^{°} \mathrm{W}$; ${27.90}^{°} \mathrm{N}, {29.96}^{°} \mathrm{N})$ with 380 cells in the east-west and 250 cells in the north-south directions, roughly. The vertical grid is divided into 10 $\sigma$-layers (11 levels) with uniform thickness. The horizontal resolution ranges between 63 m and 5 km in the cross-shore direction and between 23 m and 7.75 km in the alongshore direction. Considering that the channel width ranges between 0.25 km and 2 km, there are at least 10 grids across the narrowest channel and ~90 grids across the widest channel. The resolution progressively decreases towards offshore and away from the center of the semi-circular region (Figure 3). The southern boundary is extended offshore to the shelf edge. The fan-shaped mesh allows higher-resolution cells to be placed along the inner WLD and the zones surrounding the river mouths while reducing the grid resolution offshore and near the open boundary. The water depth in this computational region ranged between 0 and 150 m. It also included land topography with a maximum height of 19.9 m above the mean sea level, allowing the model’s wet/dry capability.

The open boundary in this model is roughly an arc around the mouth of the Atchafalaya Bay (Figure 3). The model is forced by the wind on the surface and by the tide at the open boundary using nine tidal constituents from the USACE’s Eastcoast 2001 [58] computed by ADCIRC-2DDI (i.e., the depth-integrated version of the ADCIRC; [59,60,61,62,63]).

Time-dependent and spatially uniform weather data are used as forcing on the surface. These include wind and sea-level air pressure from the station FRWL1-8766072 maintained by the National Data Buoy Center (NDBC) at Fresh Water Canal Locks, LA (Figure 1). The height of the anemometer for this station is 17.4 m. Since wind speed varies with height above the water surface and considering that the standard height used in hindcasting is 10 m, we used a logarithmic velocity profile [64] to correct the wind velocity at a different height above the water ($U_{10}$):

$$\frac{U_{10}}{U_z}={\left(\frac{10}{z}\right)}^{1/7}$$

(9)

where z is the anemometer height in m, and $U_z$ is the measured wind speed in m/s. The wind velocity is then converted to wind stress as input in the model.

The model is also provided with river discharge in the lateral boundary condition. More specifically, it includes the Wax Lake Outlet (Calumet, LA; WL, USGS station #07381590) and the Lower Atchafalaya River (Morgan City, LA; LAR, USGS station #07381590). The river discharge was unidirectional since both discharge stations are located upstream away from the estuary and hence there are no flow reversals due to tidal oscillations, although surface elevation does contain tidal signals at these stations.

The initial conditions for both surface elevation and current velocity are zero. A wet/dry scheme is included using a critical depth of $D_{min}=0.2 \mathrm{m}$. To optimize the computational efficiency, the time step is set to vary automatically based on the CFL (Courant–Friedrichs–Lewy) criterion allowing the use of a larger time step during low-speed periods, but a smaller time step when the current speed is large.

3.2. Validation: Skill Assessment

The computed water level time series at 14 stations and velocity at 4 stations are compared with the observations using the correlation coefficient (CC) and skill score (SS):

$$\mathrm{CC}=\frac{{{\displaystyle \sum }}_{i=1}^N\left(X_{mod}-\overline{X_{mod}}\right)·\left(X_{obs}-\overline{X_{obs}}\right)}{{({{\displaystyle \sum }}_{i=1}^N{\left(X_{mod}-\overline{X_{mod}}\right)}^2·{{\displaystyle \sum }}_{i=1}^N{\left(X_{obs}-\overline{X_{obs}}\right)}^2)}^{1/2}}$$

(10)

in which N is the length of the time series, $X_{mod}$ is the model output and $X_{obs}$ is the observed value.

The skill score (SS) [65] is defined as:

$$\mathrm{SS}=1-\frac{{{\displaystyle \sum }}_{i=1}^N{\left(X_{mod}-X_{obs}\right)}^2}{{{\displaystyle \sum }}_{i=1}^N{\left(X_{obs}-\overline{X_{obs}}\right)}^2}$$

(11)

Model Results with SS > 0.65 was categorized as excellent; 0.65–0.5 very good; 0.5–0.2 good; and <0.2 poor [65].

The skill assessment of the model with the combined forcing of surface wind, tide, and river discharge is calculated and summarized with the SS values shown in

Table 2. Skill score (SS) value of water level for 14 stations.

Station Number	SS Value
DL	0.5111
Big Hogs Bayou	0.6506
Small Tripod	0.6602
Large Tripod	0.7110
CSI3	0.7878
CSI6	0.7996
CSI9	0.7938
LAW	0.5603
BE	0.5227
FRWL1	0.5210
GI	0.5900
PS	0.7135
PI	0.5181
PORT FOURCHON	0.6175
WL	0.4549
LAR	0.6019

and

Table 3. Skill score (SS) value of current velocity for 4 stations (* indicates not available or applicable).

Station Number	SS Value
	Near-Surface Layer	Mid-Layer	Near-Bottom Layer
Delta #1	0.6819	*	*
Big Hogs Bayou	0.6418	*	*
CSI3	0.4155	0.6076	0.6015
CSI6	0.4032	0.4141	0.6012

. The model performance in computing water level is in the “very good” category when compared with field observations [65,66,67,68] while the model performance to compute the current velocity at the Delta #1 and Big Hogs Bayou field sampling stations is in the “excellent” category (Table 2). The correlation coefficients are generally larger than the skill scores, consistent with [67], which showed that the skill scores are generally smaller than the square of the correlation coefficients. Since we have a relatively high skill score, our correlation coefficients are even higher, and we omit the discussion for brevity.

4. Results and Discussion

4.1. Cold Front Events

Six cold front events passed the WLD region in the study period (

Table 4. Summary of cold fronts between 15 December 2012and 13 January 2013identified from meteorological observations at FRWL1 and weather maps, respectively.

Number	Time of Cold Front Passage (UTC)
	FRWL1 Observation	Weather Maps (Figures Not Shown)
1	17 December 2012, 17:48:00	17 December 2012, 21 UTC
2	20 December 2012, 11:24:00	20 December 2012, 12 UTC
3	26 December 2012, 01:30:00	26 December 2012, 00 UTC
4	28 December 2012, 21:42:00	28 December 2012, 21 UTC
5	1 January 2013, 18:48:00	1 January 2013, 21 UTC
6	10 January 2013, 11:18:00	11 January 2013, 06 UTC

; Figure 4); these events, which are associated with Migrating Cyclones (MC [3,11,13]), moved southeastward across Louisiana with average speeds of ~3 m/s. Based on the distance between Terrebonne Bay and the Atchafalaya Bay (~116 km), the average time for a cold front to pass these two locations was ~11 h.

Wind direction was variable but was dominated by southeasterly (~30%), northerly (~34%), and northwesterly (~15%) winds (Figure 4). The frequency of strong winds with speeds >10.0 m/s was rare (1.2%) and mostly associated with cold front passages coming from the northwest direction. The cold front frequency was 3 to 8 days (Table 4) and was consistent with previous studies (e.g., [1,69]). Wind patterns show that each cold front event can be identified and characterized by an abrupt change in wind direction when the surface wind quickly shifts from the prefrontal southeast wind to the stronger postfrontal northwest wind. This postfrontal wind usually persists until the next cold front approaches when the wind switches to from the southern quadrants. Each cold front event recorded during the study period lasted ~2 days except the fifth event that was characterized by a prolonged north wind lasting almost 10 days after the frontal passage. The wind speed varied significantly from ~6 m/s to ~15 m/s. Before the cold front’s arrival, air pressure decreased steadily and reached the lowest value (as low as 1000 hPa) before a rapid increase up to 1040 hPa. The air pressure usually continued to increase after the front left the area (Figure 4b) as the center of another high-pressure system moved in. At the same time, the air temperature reached the lowest value after the frontal passage (Figure 4c).

4.2. Water Flow and Transport

The WLD has seven major channels or passes radiating outward (Figure 2). At the current stage of accretion, the delta is characterized by a network of secondary channels connecting primary channels to the wetlands established at high elevations along the levees and inside the lobes/islands [24]. The structure of the channel network and connection to the inundated island interiors primarily control the water circulation and flux in this distributary system [33]. Vegetation, water level, and external forcings also modulate transport patterns [36,37].

To determine the net water flux across the delta region with a network of channels, we selected a couple dozen cross-sections (Figure 2). Water transport (or volume transport at a given time) at a given section is computed by applying the following equation:

$$F\left(x,y,t\right)={{\displaystyle \int }}_{-H\left(x,y\right)}^{\zeta }{{\displaystyle \int }}_0^L{v}_n\left(x,y,z,t\right)dldz$$

(12)

where $v_n\left(x,y,z,t\right)$ is the horizontal velocity component perpendicular to the cross-section. The total water depth is the superposition of the local undisturbed water depth, $H\left(x,y\right)$, and surface elevation, $\zeta$, from the model output. L is the width of the cross-section. The accumulated total water volume over a time period is calculated by integrating Equation (12) to yield:

$$\tilde{F}\left(x,y,t\right)={{\displaystyle \int }}_{t_0}^t{{\displaystyle \int }}_{-H\left(x,y\right)}^{\zeta }{{\displaystyle \int }}_0^L{v}_n\left(x,y,z,t\right)dldzd\tau$$

(13)

where $\tau \in \left[t_0, t\right]$ is the time interval, $t_0$ is the initial time, and final time is $t$. The subtidal water flux time series is obtained by applying the 6th-order Butterworth low-pass filter with a 0.6 cycle per day (CPD) cutoff frequency [70].

The percentage of water volume flushed out of the system during the six cold fronts ranged between 32% and 76% of the WLD’s total water volume (

Table 5. The duration, total volume transport, and percentage of volume flushed out during six cold fronts. MC: migrating cyclone.

Cold Front Event	Type	Wax Lake Delta
		Duration (h)	Percentage (%)
1	MC	41.5	32.66
2	MC	69	68.61
3	MC	44.5	76.64
4	MC	65	55.18
5	MC	65	55.18
6	MC	185	39.00

), this gives an average of 54.5% ± 16.8%. These values are larger than previously reported (Feng and Li, 2009), possibly because of the extremely shallow water in most of the WLD. It is interesting to note that the magnitude of water flushed out of the system is not necessarily dependent on the duration of the frontal events (41–185 h; Table 5). The maximum flushing out occurred during the third cold front event on 26 December 2012. This cold front had the lowest air pressure (Figure 4), maximum volume fluctuations, and maximum water elevation fluctuations (Figure 5).

This high percentage of water volume flushed out of the system underscores the significant impact of cold fronts in draining the system as water moves from the WLD to the Atchafalaya Bay before exiting to the adjacent continental shelf within a relatively short time [13]. This cold front-induced drainage is further facilitated by the extensive and shallow area where the average water column depth is ~2 m. The maximum percentage of water volume flushed out of the system reported by Feng and Li [13] by analyzing 29 cold front events was about 45%.

The simulated time-series water flow (and volume transport) in all seven channels show a diurnal variation and subtidal variations (e.g., Figure 5a) which are driven by the cold front wind. The cold front induced-water flow in the channels is significantly variable as the wind varies. Some general features, however, can be discerned as shown by the small portion of the total flow that is delivered to the Atchafalaya Bay via channelized discharge. For instance, of the accumulated water flux at these locations, only 50% of the volume at T14 flows into channels 4, 5, and 6 through T12. Further, 15% of the total flow is diverted into channels 2 and 3 through T15 and T18 while 10% moves to channel 1 via T23 (Figure 2). The remaining 25% is “lost” through the shallow water vegetated regions located among the main channels.

Based on the accumulated volumetric water flux results, all seven main channels contribute to water transport into the Atchafalaya Bay; although each channel (Figure 2b) contributes differently to the total flux. For example, the net transport from the Wax Lake Outlet transits through the cross-section T2 and then T3, which is the last channel located at the northwestern boundary. Once the water flow crosses into T4, it moves through the widest cross-section of that channel (~2 km) and then to T22 along channel 1. Channels 3, 4, 5, and 6 have the largest flow partition that is similar in range to previous flux estimates (e.g., [32]). The composite water flux channelized through channels 1 to 6 is the largest water flow occurring in the southern part of the delta. This area is also characterized by large extension of forested [71] and freshwater wetlands [29]. Thus, this area forms a complex spatial arrangement of wetland habitats flanked by subaqueous levees over which flow is driven by pressure gradients force [33].

The water flux balance can be expressed as:

$$A+B+C+D+E+F+G+Weltand Flow=TOTAL FLUX$$

(14)

where $A$, $B$, $C$, $D$, $E$, $F$, and $G$ represent the volume flux through seven main channels. The model simulations suggest that water flows into the wetland habitats account for ~25% of the total flux through the system, similar to the previous estimates by [31,72].

This water flow pattern can also be demonstrated by the water current maps showing both along- and cross-channel flows (Figure 6). The pre-frontal winds push the water up against the coastline, increasing the water level (Figure 6a,c); while the post-frontal winds push the water level down (Figure 6b,d). Flooding and water intrusion are evident in the very shallow or (sometimes) dry areas next to the channel bifurcations. Both wind magnitude and direction play a critical role in water transport across a wide topographic range. Our model simulations show significant water transport in the northwest to the southeast direction across the channels during strong post-frontal winds (Figure 6). As previously mentioned, the net water volume transported during these events had a major impact on vegetated and unvegetated areas inside the islands particularly with changes in the sediment deposition (Elliton et al. 2020), and thus needs to be considered when assessing ecological and biogeochemical processes in the delta lobes, e.g., denitrification [71] and wetland primary productivity [73]) in addition to predictions of the morphodynamic development that typically lack wind forcing, as suggested by [74].

4.3. Water Level Amplitude Spectrum

The spectrum of water levels in all field sampling stations showed a dominant diurnal tidal signal (Figure 7). Two peaks indicated the highest frequencies with values of 0.9418 CPD (or 25.48 h for one period) and 1.0009 CPD (23.79 h), close to the tidal constituents O1 (0.9295 CPD), M1 (0.9664 CPD), K1 (1.0027 CPD), S1 (1.0000 CPD), and P1 (0.9973 CPD). Although semidiurnal tidal signals were detected for most stations, their energies (see below) are lower than those observed for diurnal tides. Three major peaks were identified with frequencies of 1.884 CPD (12.74 h), 1.951 CPD (12.30 h), and 2.018 CPD (11.89 h), close to N2 (1.8960 CPD), M2 (1.9323 CPD), and S2 (2.0000 CPD). In contrast to other stations, the stations DL, BH, ST, LT, FRWL, LAW, and CSI3 (Figure 1) have relatively higher semidiurnal signals. The amplitudes of the highest semidiurnal peaks accounted for 27% to 39% of the diurnal peaks. Stations CSI3 and FRWL showed the largest semidiurnal tidal signals among all the stations. These results are consistent with other field-based water current measurements showing that areas southwest of the Atchafalaya Bay have the largest M2 tidal current across the Louisiana-Texas shelf [69]. The higher diurnal tidal components are anticipated as the region is dominated by diurnal tides [75].

The water level amplitude spectra also show periodicities longer than 2 days and are consistent with the cold front intervals. Indeed, this pattern coincides with the frequencies of cold-front-induced oscillations that are usually <0.5 CPD, and unlike tidal oscillations, cold-front-induced oscillations do not have fixed frequencies. In general, the energy is distributed throughout the low frequency range [10]. The amplitude spectra for the inshore and nearshore stations FRWL, DL, ST, LT, CSI3, BH, and LAW (Figure 1) show a strong response to wind forcing as indicated by the high-water level FFT magnitude in the 3–7 day band. In contrast, the offshore stations, e.g., CSI6, are less influenced in this range. The DL and BH spectra have similar amplitudes but slightly different patterns in the low frequency range (0–0.1 CPD) as was the case for the stations BH and LT. Further, the BH and LAW spectra have similar patterns but different amplitudes. These results indicate that although the sampling stations are in close vicinity, they are different in terms of their environment. For example, the station DL is situated in one of the main channels in the WLD and surrounded by shallow water; although BH is also in shallow waters, it is influenced by water circulation from both deltas. In contrast, the LAW station is in deeper waters close to the inner shelf (Figure 1c). There are also greater differences in velocity than in water level in response to cold front wind forcing. The greater variability in flow velocity is anticipated as the flow velocity is influenced more by the water depth.

The model results also show that the timing of the cold front is consistent with the major set-up and set-down of subtidal water level as anticipated (e.g., Atchafalaya/Vermillion Bay [13,70]; Port Fourchon [9]; and Barataria Bay [10]) (Figure 5). The variations among results from the model results (both raw results and low-pass filtered subtidal water levels) are evident in stations DL, BH, ST, LT, FRWL, LAW, CSI3, and CSI6 (Figure 5). In all cases, the pre-frontal winds from the southern quadrants caused the water level to build up. After cold fronts pass these stations, northerly winds dominate, and the water level is drawn down by the reversed wind and seaward transport of water. The variability in water level response between different events is large, ranging from <0.1 m to close to 1 m, depending on several factors: the strength and approach angle of the cold front events, wind stress, and wind direction/duration. Overall, stations DL, BH, ST, LT, FRWL, LAW, and CSI3 have relatively large variabilities compared to the other stations, apparently as a result of their relatively shallower water depths (1–5 m). The modeling results showed discernible wind-induced water level variations at all stations inside the Atchafalaya Bay where the observed water level setup was greater than 0.5 m (i.e., at DL, ST, LT, BH, FRWL, and LAW) [76] with the highest water level recorded on 26 December 2012 (Figure 5). In general, the water level variation between two successive cold front passages is caused by strong, persistent (2–3 days), and long-fetch onshore winds (>10 m/s) coming from the south and southeast quadrants. As discussed in Li et al. (2018b), the subtidal water level variation caused by cold front passages typically goes up and down over time, resembling a figure similar to a “V” or “M”.

4.4. Energy Distribution and Dominant Forcing

To evaluate the relative importance of various components of the hydrodynamics in the WLD region based on the numerical model results, we partitioned the energy into diurnal tidal, semi-diurnal tidal, subtidal, inter-tidal, and over-tidal components [10]. This is done based on the power spectrum of the water level or water velocity. The partitioning is possible because each of the five components has spectrum regimes that are mostly non-overlapping, given that our interest in the cold front effect is limited to the subtidal variations (not a short time scale variation due to, e.g., a wind gust over a few minutes or a few hours).

The results show that generally speaking, the energy partitioning varies among stations (Figure 8) but with some general patterns: the diurnal energy is generally comparable to the subtidal energy outside of the WLD region. This however is complicated by the variabilities of all components. For example, the diurnal tidal energy contribution ranges from 24–68% of the total energy while the subtidal and semidiurnal ranges from 10–46% and 10% and 21%, respectively. The diurnal tide is dominant in six out of the twelve stations followed by the subtidal component (mainly due to wind and nonlinearity). Those stations are located east of the Atchafalaya-Wax Lake complex. In contrast, the stations CSI3 and FRWL located in the western region have similar values among diurnal, semidiurnal tide, and subtidal components. What is more striking is that the energy allocation in stations within the WLD shows that the subtidal component accounted for 46–66% of total energy, which is larger than the diurnal tidal contribution of only 13–25% of the total energy. This energy partition confirms that seasonal pulsing in weather patterns—caused by the impact of atmospheric cold front passages across the Wax Lake deltaic region—is critical in regulating water circulation and the total net water flux at local scales between the delta and the Atchafalaya/Inner coastal shelf complex [10].

5. Conclusions

Results from the 3-D ECOM-si model show that both tidal and subtidal variations of water level and current velocity in the WLD and adjacent deeper open water areas are consistent with observations during a period with six cold fronts occurring between 15 December 2012 and 20 January 2013. The water transport through the delta distributary channel network under cold front weather is not solely confined within the channels; water can also have a significant lateral exchange thus hydrologically connecting the channels through the extensive wetland habitats inside the delta lobes/islands. This inter-channel lateral flow can account for approximately 25% of the total flow through the main WLD channels. During these events, the subtidal water level increased during the prefrontal period caused by the onshore southerly wind, reaching a maximum value (~1 m), and then started to drop 6–9 h prior to the frontal passage. Similarly, the subtidal influx velocity reached the highest value (~0.6 m/s) before the outward flow was developed by wind forcing from different directions after the frontal passage. The subtidal outflow reached its maximum five to eleven hours after the frontal passage because of the persistent and relatively strong winds (which can exceed 10 m/s) from the northern quadrants; the outflow, however, started to decrease as the wind weakened. The northerly winds promoted significant water flushing out of the bays as previously reported [13,70]. These cold-front-induced flushing events can last from ~40 to185 h in the WLD. Indeed, 32 to 76% of the total water in the WLD was flushed out during these cold front events analyzed in our study. This is a larger range than previously reported [13].

Finally, a spectral energy analysis revealed that both tidal and subtidal energies vary substantially. The tidal energy ranges between 20% and 70% of the total energy, while the subtidal energy varies between 10% and 45% of the total energy (Figure 8). What is most important is that within the WLD region, the weather-induced subtidal energy (46–66% of the total) is much greater than the diurnal tidal energy (13–25% of the total). This is apparently due to the shallow water in the WDL region. The shallow water (~1–2 m) implies a high nonlinear parameter (water level change divided by mean water depth) and thus a nonlinear process. The wind associated with cold fronts in winter is the main factor controlling water circulation in the WLD and is a major driver in the spatial configuration of the channel network and delta progradation rates. During the study period, river discharge had little impact on the water level and water exchange when compared to the tidal and wind effects. These results indicate that water circulation within the WLD shallow system (<2 m) is primarily dominated by wind, especially during cold fronts. Although hurricanes can cause more significant acute impacts, cold fronts are more frequent, and their integrated effect cannot be ignored when assessing hydrodynamic patterns at the regional scale; especially since cold fronts have a significant effect on shaping the WLD morphology. As discussed in [2], there were more than 1600 cold frontal events influencing the study area over a period of 40 years. The subtidal energy is more important in the WLD than in the offshore areas, thus underlying the functional role of cold fronts in this coastal region since they become more important than tides in flushing this young delta as it continues to develop and expand—even as sea level rise continues to increase in the NGOM [77,78].