Hydrodynamic and sediment transport model of the mouth of the Columbia River, Washington and Oregon, 2020-2021

A three-dimensional hydrodynamic and sediment transport model of the mouth of the Columbia River (MCR) was constructed using the Delft3D4 (D3D) modeling suite (Deltares, 2021) to simulate water levels, flow, waves, and sediment transport for the time period of September 22, 2020 to March 10, 2021. The long time period and high spatial resolution required for this study necessitated a nested modeling scheme to reduce computational expense. The nested detail model domain consists of a structured, orthogonal, curvilinear, grid that covers an area of 37 km along shore and 33 km cross shore centered on the MCR. The grid dimensions are 140 by 246 cells with resolution varying between 26 m and 1.2 km. Ten equally spaced vertical sigma layers were used to simulate 3D effects within the model domain. The grid was aligned with coastal engineering structures, including the three primary stone jetties, as well as several training dikes along the north side of the navigation channel. Flow through the structures was limited in the model by using thin dams (no transmission) or dry points. The open boundaries in the detailed model were prescribed as a time-series of Riemann invariants and salinity values for each vertical layer were derived from an overall model domain that extended roughly 150 and 100 km to the north and south of the MCR, respectively. The bathymetry for the overall, detailed, and wave grids was derived from recent data sets collected by the USGS, NOAA, and USACE between 2004 and 2020 as described in Stevens and others (2020). The source bathymetric data were converted to a common horizontal datum (NAD83) and to the land-based North American Vertical Datum of 1988 (NAVD88), then projected into the Cartesian UTM Zone 10 coordinate system (meters). Deep areas associated with the Astoria submarine canyon were removed from the model bathymetry to improve stability along the oceanic boundaries. Oceanic boundaries of the overall model were forced using astronomic tidal constituents derived from the TPXO 7.2 global tide model (Egbert and Erofeeva, 2002). A vertical offset of 1.15 m (positive values are up) derived from NOAA VDatum (version 3.2; https://vdatum.noaa.gov/) was applied at the oceanic boundary to account for the difference between local mean sea level and NAVD88. Oceanic subtidal variations were imposed at the oceanic open boundary as a time-varying correction to the astronomical tides. The subtidal time-series was derived from observations of water levels at NOAA stations 9440910 (Toke Point, WA), 9437540 (Garibaldi, OR), and 9440581 (Cape Disappointment, WA). Water-level time-series from the three stations were low-pass-filtered using a 66-hr cutoff to remove fluctuations at tidal frequencies. The low-pass-filtered values from the stations were highly correlated, and an average was applied to the oceanic model boundaries. The landward boundary of the overall model was forced with a time-series of river discharge measured at 30-minute intervals at USGS gauge 14246900 (http://waterdata.usgs.gov/usa/nwis/uv?site_no=14246900). The overall model’s oceanic and fluvial boundaries were prescribed constant salinity values of 33 and 0 practical salinity units (psu), respectively. The effects of temperature variations on circulation were neglected in the present model application.

The spectral wave model SWAN (version 41.31) was applied to both overall and detailed models to simulate waves from the continental shelf to the coastline. Wave energy was discretized into 48 frequency bins between 0.03 and 1 Hz and 36 directional bins that covered a 180-degree sector from south to north. The seaward open boundary approximately intersects the location of Coastal Data Information Program (CDIP) buoy 179 at a water depth of 181 m. The 2D, spatially uniform, time-varying energy spectra derived from measurements at CDIP buoy 179 were used to force the wave model. Space- and time-varying wind fields derived from the High Resolution Rapid Refresh (HRRR) atmospheric model (https://rapidrefresh.noaa.gov/hrrr/) were applied to simulate wind growth within the model domain. Physics in the wave model were based on ST6 (Rogers and others, 2012) and activated with a custom run script (swan.bat) that uses the find and replace text (fart.exe) utility to implement options within SWAN that are not supported by the standard Delft3d-Wave GUI. The JONSWAP bottom friction model with a coefficient of 0.038 m2s-3 and default settings for depth induced breaking and triads were included. Convergence criteria were set to 99 percent of cells and a maximum of 30 iterations to obtain full convergence for all wave cases. The overall and detailed coupled wave and flow models were run with a computational time step of 6 s to fulfill Courant stability criteria and ensure stable and accurate results. Two-way coupling between the wave model and flow model involved a nonstationary hydrodynamic calculation in combination with regular stationary wave simulations. SWAN was activated every 30 min during the hydrodynamic simulation and performed a stationary wave simulation using the water levels and depth-averaged currents passed from the flow model.

The online morphology addition to Delft3D was used to simulate sediment transport in the detailed domain at each computational time step (Deltares, 2021). Two sediment fractions representing native and berm sediment consisting of fine to medium sand were included were included. Bed stratigraphy was activated, and a total thickness of 20 m of sediment was made available throughout the model domain. The initial distributions of berm and native sediment were determined based on observed nearshore berm extent and thickness interpolated onto the computational grid. Bed level updating during the simulations was deactivated and the total mass of each sediment fraction was computed for each grid cell and vertical layer in the bed throughout the model domain.

Boundary conditions for the nested model were generated from outputs from the overall model simulation using DelftDashboard software, available as part of Open Earth Tools (https://www.openearth.nl/). The nested model simulation was broken into 4 roughly equal time periods between September 22, 2020, and March 10, 2021. This segmentation was required to reduce the amount of computer resources, specifically the amount of computer RAM, required for the simulation. After segmentation, the simulations require about 100 GB of RAM available. The time series boundary conditions (.bct) and transport conditions (.bcc) for the four segments are mcr_riv_detail_oct (Sept. 22 - Nov. 1), mcr_riv_detail_nov (Nov. 1 - Dec. 14), mcr_riv_detail_dec (Dec 14. - Jan. 26), and mcr_riv_detail_feb (Jan. 26 - Mar. 10). The nested models were run sequentially, and the successive models were initialized with conditions from the end of the prior time period.

Observations of bulk wave parameters and near-bed current velocities at approximately 2.2 m above the bed collected at the NHS buoy were compared against model predictions during low energy conditions. Qualitative and quantitative comparisons between modeled and measured water levels, wind speeds, and wave heights and periods show good agreement, with total root-mean-square errors (RMSE) of 0.11 m, 1.36 m/s, 0.22 m, and 2.21 s, respectively. However, model performance for the near-bottom currents was less accurate, with a RMSE of 0.12 and 0.14 m/s, for eastward and northward velocity components, respectively. The calibrated model was compared with observations of nearshore berm morphology and was largely able to reproduce the observed changes including the timing and magnitude of volume loss within the survey area, spreading, thinning, and onshore movement of the berm centroid. See Stevens and others (2023) for additional information.