A Wet/Dry Point Treatment Method of FVCOM, Part II: Application to the Okatee/Colleton River in South Carolina

: The wet/dry point treatment method of FVCOM was applied to simulate the tide-induced ﬂooding/drying process in the estuarine–tidal-creek–saltmarsh complex of the Okatee/Colleton River Estuary, South Carolina. The simulation results were compared with observed currents at three mooring sites and ﬂooded areas observed from remote-sensing hypsometric measurements, demonstrating that FVCOM can robustly reproduce tidal and residual currents in the river and the ﬂooding process onto the intertidal saltmarsh. The simulated ﬂow ﬁeld reveals that the Oka-tee/Colleton River Estuary is characterized by multiple residual eddies. Driven by the periodic tidal forcing, this estuarine system features a chaotic water transport process. Numerous residual eddies around the barrier complex in the Colleton River likely enhance the water exchange between the Okatee/Colleton River Estuary and the outer Broad River. A sensitivity study of ﬂooding speed to the slope of the inter-tidal zone suggests that the saltmarsh bathymetry considerably inﬂuences the water elevation near low slack water but not on the maximum water coverage area at high slack water.


Introduction
The Okatee/Colleton River Estuary (OCRE), a branch of the Broad River in South Carolina (SC), is a typically tide-controlled shallow estuary characterized by extensive intertidal saltmarshes, tidal creeks, islands, and isolated barriers ( Figure 1). The mean water depth in the main water channel of this estuary varies from~15 m in the strait between Daws Island and Colleton Neck to~2-4 m in the upstream area west of Spring Island. The tidal range exceeds 3 m (https://tidespy.com/OnePlace.php?p=1132&d=20220219&n= 1644853214, accessed on 5 July 2022, and the ratio of the tidally induced surface elevation to the mean water depth is about 0.2 in the downstream area connected to the Broad River and nearly unity in the upstream extents of the OCRE. The saltmarsh zone bounded by a 2 m water elevation line is comparable with the area of the main water channel. This area is flooded approximately twice daily by the semidiurnal M 2 tide. In 2000, a cooperative study funded by the U.S. NOAA was launched with the goal of "developing scientifically sound predictive decision-making models (tools) that integrate changes in land use patterns with effects on hydrodynamics, transport processes, and ecosystem function to assist in planning for sustainable coastal land use and resource management" [1]. The project was named "the Land Use-Coastal Ecosystem Study (LU-CES)" and the OCRE was selected as a study site for interdisciplinary field measurements and modeling efforts. Water currents and elevations were measured at three mooring sites in the OCRE over 2000-2001 ( Figure 1). Three moorings at sites A, B, and C aimed to monitor the alongriver variability of the flow and water elevation from the downstream area connected to the Chechessee and Broad rivers to the upstream end of the Okatee River. A survey of vertical current profiles with vessel-towed acoustic Doppler current profilers (ADCPs) was conducted along the main channel of the Okatee River on 8 March 2001. The R/V Gannet, a 7 m long fiberglass vessel, navigated repeated transects with 15 cycles over a 13 h period, during which current profiles were measured with a vertical bin of 0.5 m and a 15 s time interval average. Locations of the ADCP profiles were determined using a Furono-36 Differential Global Positioning System (DGPS) within an uncertainty range of 6-10 m. Harmonic analysis was used to filter the major M2 tidal constituent, and the resulting residual flow field was reconstructed along the ship track in the Okatee River [2]. To examine the water exchange processes over the estuarine-intertidal saltmarsh complex, Drs. Blanton and Andrade at Skidaway Institute of Oceanography made remote-sensing hypsometric measurements in the upstream area of the Okatee River during the flood-tidal period on 9 November 2001. The measurements were made hourly, producing six images of the water coverage area in the tidal creek-saltmarsh complex ( Figure 2). The water level was interpolated to the 5 m × 5 m square-box grid, accurately mapping the expansion of the regional flooded area during the flood-tidal period.
The field measurements made during the LU-CES study constitute a comprehensive dataset that can be used to evaluate the empirically-designed wet/dry point treatment algorithm in the Finite-Volume Community Ocean Model (FVCOM). FVCOM is an unstructured grid model developed initially by [3] and has been recently improved through a 32  Water currents and elevations were measured at three mooring sites in the OCRE over 2000-2001 ( Figure 1). Three moorings at sites A, B, and C aimed to monitor the along-river variability of the flow and water elevation from the downstream area connected to the Chechessee and Broad rivers to the upstream end of the Okatee River. A survey of vertical current profiles with vessel-towed acoustic Doppler current profilers (ADCPs) was conducted along the main channel of the Okatee River on 8 March 2001. The R/V Gannet, a 7 m long fiberglass vessel, navigated repeated transects with 15 cycles over a 13 h period, during which current profiles were measured with a vertical bin of 0.5 m and a 15 s time interval average. Locations of the ADCP profiles were determined using a Furono-36 Differential Global Positioning System (DGPS) within an uncertainty range of 6-10 m. Harmonic analysis was used to filter the major M 2 tidal constituent, and the resulting residual flow field was reconstructed along the ship track in the Okatee River [2]. To examine the water exchange processes over the estuarine-intertidal saltmarsh complex, Drs. Blanton and Andrade at Skidaway Institute of Oceanography made remote-sensing hypsometric measurements in the upstream area of the Okatee River during the flood-tidal period on 9 November 2001. The measurements were made hourly, producing six images of the water coverage area in the tidal creek-saltmarsh complex ( Figure 2). The water level was interpolated to the 5 m × 5 m square-box grid, accurately mapping the expansion of the regional flooded area during the flood-tidal period.
The field measurements made during the LU-CES study constitute a comprehensive dataset that can be used to evaluate the empirically-designed wet/dry point treatment algorithm in the Finite-Volume Community Ocean Model (FVCOM). FVCOM is an unstructured grid model developed initially by [3] and has been recently improved through a joint University of Massachusetts (UMASS)-Woods Hole Oceanographic Institution (WHOI) effort [4][5][6][7]. This model was selected for the LU-CES study because of its ability to resolve irregular coastlines, inlets, and isolated barriers, simulate the tidal flooding process over the estuarine-tidal-creek-intertidal saltmarsh complex, and ensure the volume conservation of the water transport. The model also incorporates an empirically-designed wet/dry point treatment to resolve the flooding/drying processes over intertidal saltmarshes. This method was first tested for an idealized estuarine problem [8] and then applied to the Satilla River, Georgia, with a comparison to the structured grid semi-implicit version of the Estuarine-Coastal Ocean Model (ECOMs) [8,9]. The inter-model comparisons suggested that by resolving the complex flooding/drying processes over the intertidal saltmarsh, the unstructured-grid FVCOM significantly improved the reality and accuracy of the along-channel distributions of the water elevation and residual currents in the river [9]. However, due to the lack of direct measurements of the time evolution of the flooded areas, the performance of the FVCOM wet/dry point treatment method was not evaluated by the comparison with real-time flooded area data.  [4][5][6][7]. This model was selected for the LU-CES study because of its ability to resolve irregular coastlines, inlets, and isolated barriers, simulate the tidal flooding process over the estuarine-tidal-creek-intertidal saltmarsh complex, and ensure the volume conservation of the water transport. The model also incorporates an empirically-designed wet/dry point treatment to resolve the flooding/drying processes over intertidal saltmarshes. This method was first tested for an idealized estuarine problem [8] and then applied to the Satilla River, Georgia, with a comparison to the structured grid semi-implicit version of the Estuarine-Coastal Ocean Model (ECOMs) [8,9]. The inter-model comparisons suggested that by resolving the complex flooding/drying processes over the intertidal saltmarsh, the unstructured-grid FVCOM significantly improved the reality and accuracy of the along-channel distributions of the water elevation and residual currents in the river [9]. However, due to the lack of direct measurements of the time evolution of the flooded areas, the performance of the FVCOM wet/dry point treatment method was not evaluated by the comparison with real-time flooded area data. The mooring and hypsometric data provided us with a unique opportunity to comprehensively assess the FVCOM wet/dry point treatment method, especially on its capability of capturing the current variation in the river and flooding process over intertidal saltmarsh zones in the OCRE. By comparison with observed flooded extents from the remote-sensing hypsometric measurement data, we also examined the sensitivities of flooding rate and area coverage to prescribed inter-tidal zone bathymetry. In one study [9], the FVCOM-simulated residual flow field featured multiscale eddies in the main river channel. The vessel-towed ADCP measurements verified these eddies. The multiscale eddies were also detected in the OCRE. In view of chaotic theory, when a periodic forcing was applied, the water exchange in such an estuary featured a chaotic system in which the particles moved randomly [10][11][12]. A Lagrangian particle tracking experiment was conducted to ascertain the chaotic nature of the water exchange in this estuary. The mooring and hypsometric data provided us with a unique opportunity to comprehensively assess the FVCOM wet/dry point treatment method, especially on its capability of capturing the current variation in the river and flooding process over intertidal saltmarsh zones in the OCRE. By comparison with observed flooded extents from the remote-sensing hypsometric measurement data, we also examined the sensitivities of flooding rate and area coverage to prescribed inter-tidal zone bathymetry. In one study [9], the FVCOM-simulated residual flow field featured multiscale eddies in the main river channel. The vessel-towed ADCP measurements verified these eddies. The multiscale eddies were also detected in the OCRE. In view of chaotic theory, when a periodic forcing was applied, the water exchange in such an estuary featured a chaotic system in which the particles moved randomly [10][11][12]. A Lagrangian particle tracking experiment was conducted to ascertain the chaotic nature of the water exchange in this estuary.
The remaining sections of this paper are organized as follows. In Section 2, the FVCOM model and the design of numerical experiments are described. In Section 3, the modelsimulated tidal elevation, tidal currents, and residual flow are presented and compared with the field measurements. In Section 4, the model-simulated results are compared directly with water levels and flooded areas observed from hypsometric measurements, following up with an investigation of the sensitivity of flooding speed and flooded area to the saltmarsh bathymetry. In Section 5, the Lagrangian particle tracking experiments were carried out to examine the highly complex chaotic water exchange process in the estuary. Finally, the conclusion and discussions are summarized in Section 6.

The Model
FVCOM is an unstructured grid, finite-volume, 3D primitive equation, and free surface community ocean model [3][4][5][6][7]. The FVCOM version used in this study employed a σ transformation in the vertical and non-overlapping, unstructured triangular grids in the horizontal. The governing equations were mathematically closed with the Mellor and Yamada level 2.5 turbulent closure scheme for vertical eddy viscosity [13] and the Smagorinsky eddy parameterization for horizontal diffusion coefficients [14].
The dynamics resolved by FVCOM are the same as currently popular finite-difference models such as Princeton Ocean Model (POM) [15] and Regional Ocean Model System (ROMS) [16]. A distinct difference between FVCOM and POM/ROMs is that FVCOM is solved numerically by a finite-volume flux calculation with the discretization of the integral form of the governing equations using the non-overlapping, unstructured triangular grid. The unstructured triangular grid approach in FVCOM is advantageous for resolving complex coastal boundaries and bathymetry. The flux calculation accurately represents the mass conservation in individual control volumes for water properties while retaining the coding simplicity of finite-difference methods. The discretization, time integration, and coding structure of FVCOM are described in detail in the FVCOM user manual [5,7]. A brief description of the methods and grid configuration for the OCRE is provided below.
In FVCOM, the hydrostatic primitive equations are advanced in time using a modesplitting approach in which the 2D barotropic (external) mode is integrated separately from the 3D baroclinic (internal) mode. The external mode contains the vertically integrated transport equations in which the water elevation is solved explicitly using a shorter time step constrained by min (l 1 , l 2 , l 3 ) gD, where l 1 , l 2 , and l 3 are the three side lengths of the smallest size triangle, g is gravity, D is the total local water depth, and gD is the local shallow water wave speed. A second-order accurate upwind finite-difference scheme is used for flux calculation in the integral form of the advective terms [17,18]. The modified fourth-order Runge-Kutta time-stepping scheme is selected for time integration. The internal mode consists of fully 3D governing equations and is solved using a more extended time step constrained by the phase speed of the lowest mode internal wave [3]. The linkage between external and internal modes is through the water elevation, and an external-internal mode adjustment is made at each internal time step.
The flooding/drying process in FVCOM is simulated using an unstructured wet/dry point treatment technique [19]. A viscous sublayer with a thickness of D min is added to the model to avoid the occurrence of singularity when the local water depth approaches zero. In this system, the wet or dry criterion for node points is given as and for triangular cells where h B is the bathymetric height related to the river edge where the mean water depth is zero.î,ĵ, andk are integer numbers to identify the three-node points of a triangular cell. When a triangular cell is treated as dry, the velocity at the centroid of this triangle is specified to be zero, and no flux is allowed through the three side boundaries of this triangle. This triangular cell is removed from the flux calculation in the control volume.
A detailed discussion of this unstructured grid wet/dry point treatment technique is given in [19]. For brevity, it is not included in this paper.

Design of Numerical Experiments
The numerical experiment was conducted using the nested OCRE and Okatee tidal creek (OTC) FVCOMs with grid configuration shown in Figure 3. The domains were meshed using a non-overlapping unstructured triangular grid, with the horizontal resolution ranging from~250 m in the main channel of the estuary to~25 m upstream of the OCRE in the OCRE-FVCOM and from~5 m in the main river channel to~10 m over the intertidal saltmarsh in the OTC-FVCOM. Two open boundaries spanned the straits between Lemon and Daws Islands and between Daws Island and Colleton Neck in the OCRE-FVCOM. A uniform sigma coordinate transformation with 11 levels was used in the vertical, corresponding to a vertical resolution of 0.5 m or less in the computational regions. The time step was 1.242 s for the external mode for the OCRE-FVCOM and 0.2484 s for the OTC-FVCOM. In both domains, the external-internal time step ratio was 10. This ratio was determined based on the stability analysis results presented in [19]. The 5 m resolution topographic data obtained from hypsometric measurements were used to configure the OTC-FVCOM, allowing to make an Apple-to-Apple comparison between simulated and observed flooded areas during the flood-tidal period.
where ℎ is the bathymetric height related to the river edge where the mean water depth is zero. , , and are integer numbers to identify the three-node points of a triangular cell. When a triangular cell is treated as dry, the velocity at the centroid of this triangle is specified to be zero, and no flux is allowed through the three side boundaries of this triangle. This triangular cell is removed from the flux calculation in the control volume. A detailed discussion of this unstructured grid wet/dry point treatment technique is given in [19]. For brevity, it is not included in this paper.

Design of Numerical Experiments
The numerical experiment was conducted using the nested OCRE and Okatee tidal creek (OTC) FVCOMs with grid configuration shown in Figure 3. The domains were meshed using a non-overlapping unstructured triangular grid, with the horizontal resolution ranging from ~ 250 m in the main channel of the estuary to ~ 25 m upstream of the OCRE in the OCRE-FVCOM and from ~ 5 m in the main river channel to ~ 10 m over the intertidal saltmarsh in the OTC-FVCOM. Two open boundaries spanned the straits between Lemon and Daws Islands and between Daws Island and Colleton Neck in the OCRE-FVCOM. A uniform sigma coordinate transformation with 11 levels was used in the vertical, corresponding to a vertical resolution of 0.5 m or less in the computational regions. The time step was 1.242 s for the external mode for the OCRE-FVCOM and 0.2484 s for the OTC-FVCOM. In both domains, the external-internal time step ratio was 10. This ratio was determined based on the stability analysis results presented in [19]. The 5 m resolution topographic data obtained from hypsometric measurements were used to configure the OTC-FVCOM, allowing to make an Apple-to-Apple comparison between simulated and observed flooded areas during the flood-tidal period.  The OCRE-FVCOM was driven by the tidal forcing at two open boundaries connected to the Broad River Estuary. Tidal forcing at open boundaries was provided by a regional coarse-resolution estuarine model for the Broad River Estuary (the model grid was not included in this manuscript). Amplitudes and phases of five major tidal constituents (M 2 , N 2 , S 2 , O 1 , and K 1 ) at boundary nodes were specified directly using the output of the Broad River Estuary model. The tidal boundary forcing was adjusted to produce the best match between model-simulated and observed sea levels at available measurement sites. The primary freshwater discharge was from the upstream area of the Broad River. The river runoff from the upstream end of the OCRE was too minimal to be accounted for. It was ignored in the OCER-FVCOM run, while it was included in the sub-domain OTC-FVCOM experiments. The river discharge added into the OTC-FVCOM was about 0.01 m 3 /s [20]. The OTC-FVCOM shared the same common cells with the OCRE-FVCOM at the nesting boundary. The OCRE-FVCOM provided the forcing to drive the OTC-FVCOM through the nesting boundary.
Numerical experiments were carried out following three steps. The first was focused on tidal simulation, the second on the flooding/drying process, and the third on particle tracking. The model results were validated through direct comparisons with field measurements. The discussion on the chaotic water transport process was primarily based on model-computed particle trajectories. The LU-CES field measurements showed that the water in the OCRE was vertically well mixed throughout the year. For this reason, no salinity or temperature fields were included in our present numerical experiments.

The Water Elevation
The OCRE-FVCOM was run for more than 30 days (60 tidal cycles) to resolve the spring-neap tidal cycle. The M 2 , S 2 , N 2 , and O 1 , and K 1 were separated from the modelpredicted surface elevation and current data using Foreman's tidal harmonic analysis program [21]. Model-predicted amplitudes and phases of these five major tidal elevations were directly compared with the observational data at three mooring sites (numbered A, B, and C) shown in Figure 1 26 April 2001. A 3 h low-passed filter was used to re-sample these data at an hourly interval. It was the state-of-the-art method to process the tidal measurement data. In this method, resampling did not change tidal energy spectra within an independent time scale of the tidal signal [22][23][24]. The same Foreman's harmonic analysis program was used to calculate the amplitudes and phases of observational tidal constituents.
The OCRE-FVCOM provided a robust simulation of M 2 , S 2 , N 2 , K 1, and O 1 tidal elevations. The model predicted that the M 2 tidal constituent dominated the tidal motion in the OCRE. The S 2 or N 2 to M 2 ratio was about 0.2, and K 1 or O 1 to M 2 was about 0.1. The amplitudes of the five tidal elevation constituents increased gradually from the Colleton River to the Okatee River, with a net upstream-downstream difference of about 10 cm in M 2 , <5 cm in S 2 and N 2 , and <3 cm in K 1 and O 1 . The model-predicted spatial trend matches the observed data at sites A, B, and C ( Figure 4). A quantitative comparison of amplitudes and phases for M 2 , S 2 , N 2 , K 1, and O 1 tidal elevations was listed in Table 1. The error bars were estimated by the error analysis program from [25]. The overall differences between simulated and observed amplitudes and phases were 1.3 cm and 0.3 • for M 2 , 0.3 cm and 8 • for S 2 , 1.2 cm and 7 • for N 2 , 0.3 cm and 5 • for K 1 , and 0.1 cm and 5 • for O 1 , respectively.
The tidal simulation results, however, showed a relatively large bias in the tidal phase at sites B and C, especially for S 2 and N 2 tidal constituents. This phase error increased from site B to site C. To investigate a possible reason for this error, we examined the sensitivities of the computed amplitude and phase to the turbulent mixing parameterization. The results showed that the tidal phase was more sensitive to bottom roughness (z o ) specified in the model than tidal amplitude. One study [26] examined the tidal simulation uncertainty over Georges Bank due to turbulence mixing parameterization. They found that tidal mixing was generally predominated in the coastal region by a local 1-dimensional balance between turbulent shear production and dissipation. The spatial distribution of local turbulent dissipation could be influenced significantly by the spatial-bathymetric resolution and spatial distribution of z o . Since no observational data were available for the spatial distribution of z o in this estuary, a constant value of 0.001 m was specified everywhere in the model. The relatively large errors found in the model-predicted phases of S 2 and N 2 at sites B and C in the Okatee River were probably due to the z o uncertainty and bathymetry in the shallower muddy upstream region.
(zo) specified in the model than tidal amplitude. One study [26] examined the tidal simulation uncertainty over Georges Bank due to turbulence mixing parameterization. They found that tidal mixing was generally predominated in the coastal region by a local 1dimensional balance between turbulent shear production and dissipation. The spatial distribution of local turbulent dissipation could be influenced significantly by the spatialbathymetric resolution and spatial distribution of zo. Since no observational data were available for the spatial distribution of zo in this estuary, a constant value of 0.001 m was specified everywhere in the model. The relatively large errors found in the model-predicted phases of S2 and N2 at sites B and C in the Okatee River were probably due to the zo uncertainty and bathymetry in the shallower muddy upstream region.

Tidal Currents
The bottom-mounted RDI 1200 kHz Work Horse ADCPs were deployed at sites B and C during the same period as the measurement of the tidal elevation at those two sites.
A 0.5 m bin-averaging was used for the ADCP measurements. The quality data were received only within the water column away from the transducer by a distance required for the instrument design to the processing echo sound data [27]. At these two sites, the ADCP vertical velocity profiles were available in the water column about 2 m below the surface and above the bottom. The Foreman's tidal harmonic analysis program was used to calculate the tidal current ellipse profiles [21]. The model-simulated vertical profiles of M 2 , S 2 , N 2 , K 1 , and O 1 tidal currents were in good agreement with the ADCP-measured tidal currents at site B ( Figure 5) and site C ( Figure 6). A detailed comparison between model-simulated and observed major axes, minor axes, orientations, and phases is shown in Table 2. Considering the uncertainty in current measurements shown in Table 2, the model reasonably reproduced the vertical shear of the major axis of tidal currents at both sites B and C. Relatively large errors in the current phase and orientation found at the measurement sites were believed to result from the incorrect z o parameterization in the muddy regions of the estuary. As previously mentioned, the sensitivity of tidal current phases and orientation to z o was recognized for a long time in the coastal ocean community [26]. Since z o was a temporospatial function relating to the seabed types, the experiments by tuning a spatially uniform z o would not enable to improve an overall tidal current simulation over such an estuarine-tidal-creekintertidal saltmarsh complex.

Residual Currents
The residual flow was defined as the averaged current after removing tidal currents constructed by the five major tidal constituents. The tidal currents were extracted using Foreman's harmonic analysis method [21]. In the OCRE, the M 2 tidal currents dominated. The residual current patterns predicted by filtering the five major tidal constituents over a long time series showed no significant difference from the residual current computed in the model run forced solely by the M 2 tidal component.
The model results showed that residual currents in the OCRE were characterized by multiple clockwise and anticlockwise eddy-like flows (Figure 7). In the main river channel, where no islands or barriers exist, eddies were centered on the main axis of the river, with their velocities increasing towards the coast. The maximum speed of a residual vortex reached 25 cm/s near the surface and~15 cm/s near the bottom. Between sites A and B in the Colleton River, where several islands and barriers were located, the eddy shape was constrained by proximal island and barrier geometries. Around those areas, the spatial structure of the residual currents was much more complex. This finding was consistent with our previous model results conducted in the Satilla River Estuary, Georgia, where the residual flow was also predominated by multi-eddies [9]. Those eddies were verified by vessel-towed ADCP measurements [9].

Residual Currents
The residual flow was defined as the averaged current after removing tidal currents constructed by the five major tidal constituents. The tidal currents were extracted using Foreman's harmonic analysis method [21]. In the OCRE, the M2 tidal currents dominated. The residual current patterns predicted by filtering the five major tidal constituents over a long time series showed no significant difference from the residual current computed in the model run forced solely by the M2 tidal component.
The model results showed that residual currents in the OCRE were characterized by multiple clockwise and anticlockwise eddy-like flows (Figure 7). In the main river channel, where no islands or barriers exist, eddies were centered on the main axis of the river, with their velocities increasing towards the coast. The maximum speed of a residual vortex reached 25 cm/s near the surface and ~ 15 cm/s near the bottom. Between sites A and B in the Colleton River, where several islands and barriers were located, the eddy shape was constrained by proximal island and barrier geometries. Around those areas, the spatial structure of the residual currents was much more complex. This finding was consistent with our previous model results conducted in the Satilla River Estuary, Georgia, where the residual flow was also predominated by multi-eddies [9]. Those eddies were verified by vessel-towed ADCP measurements [9]. A vessel-towed ADCP survey was conducted along the Okatee River on March 1, 2001. The ADCP was mounted on a sled and towed on one side of the R/V Gannet. The measurements were made during daylight and lasted for 13 h. The vessel was kept at a steady speed at ~ 5-6 knots, and velocity was recorded with a vertical interval of 0.5 m and a time interval of 15 s. Harmonic statistical analysis extracted the amplitude and phase of the M2 tidal velocity. A vertically averaged residual velocity field was constructed after removing tidal currents, with uncertainty in a range of about 4-16 cm/s. This uncertainty was mainly caused by the measurement bias due to the variation of irregular water depths along the ship track and the finite vessel speed [2]. A vessel-towed ADCP survey was conducted along the Okatee River on 1 March 2001. The ADCP was mounted on a sled and towed on one side of the R/V Gannet. The measurements were made during daylight and lasted for 13 h. The vessel was kept at a steady speed at~5-6 knots, and velocity was recorded with a vertical interval of 0.5 m and a time interval of 15 s. Harmonic statistical analysis extracted the amplitude and phase of the M 2 tidal velocity. A vertically averaged residual velocity field was constructed after removing tidal currents, with uncertainty in a range of about 4-16 cm/s. This uncertainty was mainly caused by the measurement bias due to the variation of irregular water depths along the ship track and the finite vessel speed [2].
Despite a relatively large measurement error, the along-river distribution of the vesseltowed ADCP-derived residual current agreed qualitatively with the model residual currents in the same plane (Figure 8). For example, the observed residual currents showed a net outward flow of 5-10 cm/s along the right side of the Okatee River. This flow was captured in the distribution of the residual current predicted by the model. Since the measurement uncertainty was larger than the residual flow in the weak velocity areas, we did not quantitatively compare the simulated and observed velocities long the ADCP survey tracks.
Okatee River. Based on the ADCP measurements, it appeared that there was a net outward flux in the Okatee River, which violated basic volume conservation principles. The model resolved the flow field in the under-sampled regions, showing that a return flow on the left side of the river compensated for the water loss due to the outflow on the right side of the river. It suggested that the observed ADCP current data could be used to evaluate the model performance. Meanwhile, due to the flow complicit in this estuary, the model could also be used to interpret the observed finding.

The Flooding/Drying Process
The 3D wet/dry point treatment used in FVCOM produced the flooding/drying process in the estuarine-tidal-creek-intertidal saltmarsh complex over tidal cycles. Figure 9ad showed the distributions of the sea level and vertically averaged tidal velocity vectors at the maximum flood, the transition from flood to ebb (high water), the maximum ebb, and the ebb to flood transition (low water). Water flooded onto the OCRE from two paths around Spring Island from the Broad River. At high water, most of the intertidal saltmarsh area was covered by the water. In turn, the water drained back during the ebb period and sequestered in the main channel at low water.
One of the critical issues in this modeling application was to verify the accuracy of the simulation of the water transport during the flooding/drying process over inter-tidal zones. Many factors, such as the poor resolution of bathymetry over the saltmarsh, model horizontal resolution, turbulence mixing parameterization, model time steps, etc., could alter the accuracy of the model simulation. As previously mentioned, the remote-sensing hypsometric curves measured in the upstream tidal creek area of the Okatee River during the flood-tidal period provided unique high-resolution image sets to assess the FVCOM This comparison provides a good example of the need to integrate model and observed data. Due to limited coverage of the survey area, the vessel-towed ADCP along the cruise track was insufficient to resolve the 3D spatial structure of the residual flow in the Okatee River. Based on the ADCP measurements, it appeared that there was a net outward flux in the Okatee River, which violated basic volume conservation principles. The model resolved the flow field in the under-sampled regions, showing that a return flow on the left side of the river compensated for the water loss due to the outflow on the right side of the river. It suggested that the observed ADCP current data could be used to evaluate the model performance. Meanwhile, due to the flow complicit in this estuary, the model could also be used to interpret the observed finding.

The Flooding/Drying Process
The 3D wet/dry point treatment used in FVCOM produced the flooding/drying process in the estuarine-tidal-creek-intertidal saltmarsh complex over tidal cycles. Figure 9a-d showed the distributions of the sea level and vertically averaged tidal velocity vectors at the maximum flood, the transition from flood to ebb (high water), the maximum ebb, and the ebb to flood transition (low water). Water flooded onto the OCRE from two paths around Spring Island from the Broad River. At high water, most of the intertidal saltmarsh area was covered by the water. In turn, the water drained back during the ebb period and sequestered in the main channel at low water.
One of the critical issues in this modeling application was to verify the accuracy of the simulation of the water transport during the flooding/drying process over intertidal zones. Many factors, such as the poor resolution of bathymetry over the saltmarsh, model horizontal resolution, turbulence mixing parameterization, model time steps, etc., could alter the accuracy of the model simulation. As previously mentioned, the remotesensing hypsometric curves measured in the upstream tidal creek area of the Okatee River during the flood-tidal period provided unique high-resolution image sets to assess the FVCOM performance, especially its capability of simulating a natural flooding process in this estuary.
To validate the reality and accuracy of the FVCOM wet/dry treatment method, we first ran the OCRE-FVCOM using the real-time tidal forcing consisting of five major tidal constituents over 1-16 November 2001, and then the refined-grid OTC-FVCOM through the nesting approach over the same period. The simulated sea levels and flooded areas were directly compared with the hypsometric data over the measurement period from 9:30 a.m. to 3:04 p.m. on 9 November 2001. It should be pointed out that there was no bathymetry data over intertidal saltmarsh areas covering the entire OCRE-FVCOM domain. Over those areas, the model grid was configured using a linear function from the edge of the river's main channel to the 2 m water elevation. Despite it, the OCRE-FVCOM was robust in reproducing the evolution of the water level and flooded area in the OCRE except near the transition time from ebb to flood (Figures 10 and 11). The modelsimulated maximum flooded area at high water was very close to that observed from the hypsometric measurements.

2, 10, x FOR PEER REVIEW
14 of 20 performance, especially its capability of simulating a natural flooding process in this estuary. To validate the reality and accuracy of the FVCOM wet/dry treatment method, we first ran the OCRE-FVCOM using the real-time tidal forcing consisting of five major tidal constituents over 1-16 November 2001, and then the refined-grid OTC-FVCOM through the nesting approach over the same period. The simulated sea levels and flooded areas were directly compared with the hypsometric data over the measurement period from 9:30 a.m. to 3:04 p.m. on 9 November 2001. It should be pointed out that there was no bathymetry data over intertidal saltmarsh areas covering the entire OCRE-FVCOM domain. Over those areas, the model grid was configured using a linear function from the edge of the river's main channel to the 2 m water elevation. Despite it, the OCRE-FVCOM was robust in reproducing the evolution of the water level and flooded area in the OCRE except near the transition time from ebb to flood (Figures 10 and 11). The model-simulated maximum flooded area at high water was very close to that observed from the hypsometric measurements.    The comparison results exceeded our expectations, as there was no observed bathymetric data over the intertidal saltmarsh areas available to configure the OCRE-FVCOM. To test the sensitivity of the geometric slope over the intertidal saltmarsh zone to the water level and flooded area, we re-ran the OCRE-FVCOM with a slope specified using a parabolic function from the edge of the river's main channel to the 2 m water elevation (see Figure 12). The results indicated that the bathymetry shape over the intertidal zone did change the tidal amplitude but not the phase. The change mainly affected the speed of flooding and thus the simulated flooded area. The area was altered significantly by the shape of the slope during the early flooded period but had little influence on the maximum flooded area and maximum water level at high water. This might explain why OCRE-FVCOM provided a reasonable simulation for the maximum flooded area and water level at high water but not for the flooded area in the early flooding period. On the other hand, good comparison results shown in Figures 10 and 11 suggested that the OCRE-FVCOM was accurate enough to simulate the flooding process in the OCRE regarding the water exchange over the estuarine-tidal-creek-intertidal saltmarsh complex. The comparison results exceeded our expectations, as there was no observed bathymetric data over the intertidal saltmarsh areas available to configure the OCRE-FVCOM. To test the sensitivity of the geometric slope over the intertidal saltmarsh zone to the water level and flooded area, we re-ran the OCRE-FVCOM with a slope specified using a parabolic function from the edge of the river's main channel to the 2 m water elevation (see Figure 12). The results indicated that the bathymetry shape over the intertidal zone did change the tidal amplitude but not the phase. The change mainly affected the speed of flooding and thus the simulated flooded area. The area was altered significantly by the shape of the slope during the early flooded period but had little influence on the maximum flooded area and maximum water level at high water. This might explain why OCRE-FVCOM provided a reasonable simulation for the maximum flooded area and water level at high water but not for the flooded area in the early flooding period. On the other hand, good comparison results shown in Figures 10 and 11 suggested that the OCRE-FVCOM was accurate enough to simulate the flooding process in the OCRE regarding the water exchange over the estuarine-tidal-creek-intertidal saltmarsh complex. The sub-domain OTC-FVCOM was configured with the high-resolution bathymetry obtained from hypsometric measurements. The 5 m resolution bathymetric data used in this model covered the Okatee tidal creek and intertidal saltmarsh. The OTC-FVCOM was capable of reproducing the expanding process of the observed flooded areas during the flood-tidal period (Figure 13), robustly matching with the observed flooded area displayed in Figure 2. Both the observation and model consistently showed that the flood tidal period lasted about 6 h. The water over the intertidal saltmarsh was almost com- The sub-domain OTC-FVCOM was configured with the high-resolution bathymetry obtained from hypsometric measurements. The 5 m resolution bathymetric data used in this model covered the Okatee tidal creek and intertidal saltmarsh. The OTC-FVCOM was capable of reproducing the expanding process of the observed flooded areas during the flood-tidal period (Figure 13), robustly matching with the observed flooded area displayed in Figure 2. Both the observation and model consistently showed that the flood tidal period lasted about 6 h. The water over the intertidal saltmarsh was almost completely drained back to the Okatee River at the lowest water level (9:38 a.m.) (Figure 13a). The flood tide gradually advected the water into tidal creeks during the subsequent 1.5 h, and significant flooding onto the intertidal saltmarsh occurred within 2 h after that. The agreement between observed and simulated flooded areas at six measurement times suggested that for the given accurate bathymetry, the FVCOM wet/dry point treatment method was robust to simulate the flooding process over the estuarine-tidal-creek-intertidal saltmarsh complex. This finding was also demonstrated in the FVCOM applications to the Plum Island Sound, MA [27], storm-induced coastal inundation in Scituate Harbor and Massachusetts Bay [28][29][30], and tsunami-induced coastal inundation over the Japanese coast [31].

Chaotic Nature of Water Transport
The model-simulated residual flow field showed that the OCRE was charact

Chaotic Nature of Water Transport
The model-simulated residual flow field showed that the OCRE was characterized by multiple residual eddies. These residual eddies usually acted as retention zones or dynamic barriers to trap biological and chemical tracers. When an extra periodic tidal forcing was applied, however, the particle trajectories became chaotic, and rapid water exchanges occurred between the OCRE and Broad River. The feature was consistent with the chaotic theory described in [10][11][12].
By releasing 100 neutral buoyancy fluid particles throughout the water column in the upstream area of the Okatee River, we traced the trajectories of particles and the variance change of particles as a mass in the estuary. The purpose of this experiment was to apply the Fickian model to estimate the horizontal diffusion coefficient based on the change of variance of the mass consisting of particles [32]. The experiments failed due to random movements of particles. An example of randomly-moving particle trajectories is shown in Figure 14. The particle released in the upstream area of the Okatee River moved outward during the first eight tidal cycles. When it arrived in the Colleton River, where strong closed residual eddies were located around the barrier complex, it suddenly moved outward faster than the advective tidal current speed. When we released particles around those barrier complexes in the Colleton River, trajectories of particles become chaotic in the first tidal cycle. These features were also observed in the trajectories of particles released near the bottom. This result suggested that the water exchange process in the OCRE was chaotic in nature. The existence of barrier complexes acted as an accelerator to enhance the water exchange between the OCRE and the outer Broad River. This chaotic system played an essential role in maintaining a healthy environment in this estuary. was chaotic in nature. The existence of barrier complexes acted as an accelerator to enhance the water exchange between the OCRE and the outer Broad River. This chaotic system played an essential role in maintaining a healthy environment in this estuary. The particle tracking experiment also raised a fundamental issue on kinematics related to the chaotic process in estuaries. The model results showed that multiple residual eddies also characterized the residual flow in the upstream area of the Okatee River, but under the same periodic tidal forcing environment, no chaotic features of particle movements were found in that area. It suggested that the existence of double circulation cells and periodic tidal forcing were necessary but not sufficient conditions to cause the chaotic features. The fact that the chaotic process occurred in the river segment with a broader width around the barrier complex in the Colleton River implied that this process might The particle tracking experiment also raised a fundamental issue on kinematics related to the chaotic process in estuaries. The model results showed that multiple residual eddies also characterized the residual flow in the upstream area of the Okatee River, but under the same periodic tidal forcing environment, no chaotic features of particle movements were found in that area. It suggested that the existence of double circulation cells and periodic tidal forcing were necessary but not sufficient conditions to cause the chaotic features. The fact that the chaotic process occurred in the river segment with a broader width around the barrier complex in the Colleton River implied that this process might be related to the geometry and bathymetry of the OCRE.

Summary and Discussion
The tide-induced flooding/drying process over the estuarine-tidal-creek-intertidal saltmarsh complex of the OCRE, SC, was simulated by FVCOM. The model results clearly showed that the OCRE was characterized by the relatively strong M 2 tidal motion overlapped with multiple residual eddies. Driven by periodic tidal forcing, the water exchange process in the estuary was chaotic in nature. These numerous residual eddies in the Colleton River could enhance the water exchange between the OCRE and the outer Broad River.
The data-model comparisons demonstrated that FVCOM could reproduce the tidal motion, residual flow field, and the water flooding onto the intertidal saltmarsh in the OCRE. A sensitivity study suggested that the details of the bathymetry over the intertidal saltmarsh could considerably influence the water elevation and the flooded area near the ebb to flood tidal transition, but not on the maximum water coverage area at the high slack water. For given accurate bathymetry, the FVCOM wet/dry point treatment method was robust to simulate the flooding processes onto the intertidal saltmarsh.
It should be noted that the amplitude of the OCRE-FVCOM-simulated M 4 tidal constituent, a sign of the strength of nonlinear tidal interaction, was much weaker than the observed value estimated from tidal measurements at mooring sites. The modeling experiments using the OTC-FVCOM suggested that a strong M 4 tidal current was a crucial factor in the asymmetric feature of the tidal motion over the intertidal saltmarsh area in the tidal creek [20]. The failure to reproduce the strength of the M 4 tidal constituent in the OCRE-FVCOM experiment was believed to be a result of poor bathymetry resolution over the intertidal saltmarsh area and insufficient model resolution in the main estuary. This issue was addressed in Chen et al. [20].