Coupling a Parametric Wave Solver into a Hydrodynamic Circulation Model to Improve Efﬁciency of Nested Estuarine Storm Surge Predictions

: Efﬁciency in storm surge modeling is crucial for forecasting coastal hazards in real-time. While computation cost may not be the main concern for organizations with ample resources, the robustness of forecasts generated by most parties are restricted by wall-clock time. The Parametric Wave Solver (PARAM) was developed by Boyd and Weaver (2021) as an alternative to computationally expensive wind–wave models when modeling restricted estuarine environments. For this study, PARAM has been tightly coupled with the ADCIRC hydrodynamic model to create ADCparam, then integrated into the Multistage mini-ensemble modeling system (MMEMS), a one-way nesting framework for modeling waves and circulation in coastal estuaries developed by Taeb and Weaver (2019). In the MMEMS framework, ADCIRC + SWAN is used to simulate the coarser ocean domain and ADCparam is applied to the nested high resolution estuarine mesh. ADCparam has greatly reduced computation time for the high resolution nested sub-model compared to the third-generation wave model originally used. While the PARAM wave solution shows dissimilarities with the SWAN solution, signiﬁcant wave height and wave period results are consistent and warrant further pursuit of the parametric wave ensemble method as a substitute to SWAN within MMEMS. ADCparam models demonstrated run times up to 51% faster than ADCIRC coupled with SWAN, an established iterative wave model tightly coupled to ADCIRC and packaged with the MMEMS repository. ADCparam wall time is comparable to running ADCIRC without wave forcing, nearly eliminating the computational cost of including the wave forcing in the high-resolution estuarine domain of MMEMS. Computational efﬁciency is greatly increased while maintaining solution integrity. Though ADCparam, and its application to MMEMS, are still being reﬁned and validated, the coupled model system has proven to be an efﬁcient, viable path for implementing waves in any estuarine circulation model.


Introduction
Oceanographers and engineers use numerical modeling to understand physical processes related to coastal morphology, hydrodynamic flow, waves, and turbulence [1]. Although real-world observation is necessary to derive descriptions of these phenomena and their interactions, numerical computer models extend the range of insights. Using numerical schemes converted from differential relationships [2], computers are capable of calculating simulations for entire oceans. By simulating the hydrodynamic response of the ocean to known (hindcasted) or predicted (forecasted) tidal and atmospheric conditions, researchers can recreate or make predictions of water levels, wave spectra, currents, Lagrangian transport, turbulence, and/or bathymetric morphology at every point in a basin-scale domain [1]. Applications for these models include (but are not limited to) coastal flooding, maritime safety, beach nourishment, dredging, marine spatial planning, and climatology. Protecting lives and infrastructure near coastlines creates the necessity for generating accurate predictions as efficiently as possible. This need is emphasized by the increasing magnitude of severe storms [3]. This study seeks to improve real-time storm surge hazard forecasts for estuarine environments with limited computational power. By expediting computation of wind-driven waves and their contribution to the total forcing in a hydrodynamic model, overall computational wall-clock time can be significantly reduced.

Wind-Wave Hydrodynamics
Ocean waves are categorized by their height, period, and the water depth in which they are propagating. Wave generating mechanisms (e.g., wind, tide, earthquakes, etc.) are associated with creating waves in a characteristic frequency range. Wind-generated waves are the primary type pertaining to this study, and have periods ranging from fractions of a second, up to about five minutes for infragravity waves in the open ocean [4] (Figure 1). Most numerical modelers are concerned with waves with periods greater than one second, and less than 30 s [5]. Waves in this range of periods are called wind-generated gravity waves, as gravity is the restoring force pulling the wave crests downward. Surface tension takes over as the primary restoring force for waves with periods much lower than one second, and are regarded as capillary waves [6]. Long period waves, with periods ranging from minutes to days, Figure 1, are also included in the numerical models as periodic forcing frequencies, which mainly simulate tides [5]. Wave-solving programs are generally coupled to hydrodynamic models in order to include the effects of waves generated by the incident wind field on circulation [9]. While gravity is the restoring force for these waves, the initial disturbance due to wind is not as clearly defined. As wind blows over a calm ocean, there is friction between the two mediums. Waves form from the variation in pressure forces and turbulent processes, which amplify as the water surface grows more uneven [4]. Within the wave generation zone, wave heights and periods are dependent on the wind speed, the duration that the wind blows, and the distance over which the wind is blowing (i.e., fetch).
Since waves are not linear in nature, they drive currents and mass transport through turbulent processes and variations in the local pressure field. As waves propagate, shoal, and break unevenly over complex bathymetry, water particles follow the ensuing pressure gradients. Wave breaking alone raises water levels through set-up, caused by the momentum flux between the breaking wave and the water. The kinetic energy from the dissipating waves is converted to potential energy in the water column, raising the water elevation, as well as kinetic energy in the form of turbulence and longshore currents [10]. The increase of mean water level at the coast can account for as much as 30% of the total water level rise, or storm surge [9]. The inclusion of wave physics in circulation models is important for accurate calculation of both water elevation and currents.

Modeling Solutions
Due to computational limitations, the governing equations of fluid mechanics cannot be solved analytically for realistic sized domains, so they are discretized with approximate algebraic expressions. These expressions are solved at a finite number of domain nodes or elements for a fixed rate of time [2]. This method of computational fluid dynamics (CFD) has become widely applied since computers powerful enough to solve them have become available, and has allowed for hydrodynamic processes to be resolved on global geographic scales. Verified systems of CFD programs, such as ADCIRC [11], Delft3D [12], SLOSH [13], CMSFLOW [14], FVCOM [15], ROMS [16], and CH3D-SMSS [17], are a few of the most notable packages for modeling basin-scale ocean processes. When inputs, computational techniques, and parameterizations are equal between models, there is little difference in solutions between the solutions computed by the accepted circulation models [18].
Hydrodynamic modeling programs are fundamentally distinguished by their numerical schemes and physical parameters. ADCIRC (Advanced CIRCulation), the program package selected as a foundation for this study, applies a finite element spatial scheme to a modified form of the shallow water wave continuity equation [11].
ADCIRC starts with a flexible mesh grid space of triangular elements, each with three nodes storing coordinates, elevation, and frictional characteristics [11]. An unstructured mesh can conform to bathymetry for smoother gradients between elements, ensuring continuity and stability across the model [19]. Tides in the form of harmonic constituents are forced at open boundaries to initiate the finite element sweep, typically with a ramp period to gradually increase amplitudes, while 10 m wind and atmospheric pressure fields are applied to corresponding nodes for each timestep.
Similar numerical processes are used to compute wave spectra across ocean domains, which can be coupled with hydrodynamic models to account for wave action contribution to the overall water level. Radiation stress gradients from wave breaking can be an order of magnitude greater than gradients from wind stress alone [20]. SWAN (Simulating WAves Nearshore), one of the wave models routinely coupled with ADCIRC, solves for wave energy density using an implicit discretization of the modified action balance equation, categorizing it as a third-generation wave model [21,22]. Phase averaging third-generation wave models brought improvements to wave solutions under extreme wind fields when they were first developed in 1988 [23]. Featured advancements include Komen's energy dissipation source function [24], and Hasselmann's discrete interaction [25]. SWAN is set apart from other third-generation models, such as STWAVE and WAVEWATCH III, by the inclusion of processes unique to shallow water environments, making it applicable to shallow restricted estuaries [26][27][28][29].
Although the third-generation wave model is markedly stable, the drawback to this type of scheme is prolonged computational time. In this formula, wave energy is a function of five dimensions: time, two axes of geographic space, and two axes of spectral space [21,22]. As such, each time SWAN is called in the coupled model (i.e., every 600 timesteps), it must iterate for values of wave parameters that are dependent on multiple preceding values in space and time. Due to the fact that the solution can take considerable time to fully converge, SWAN input includes the option to limit the number of iterations.

Background
A critical requirement of real-time coastal hazards modeling is calculating results in a timely manner. Forecasters and emergency managers need model solutions early enough to incorporate the results in their analysis and provide ample warning to those at risk. Computation time can drive the utility of the guidance product. For engineering consultant firms and local government agencies limited by computer processing power, one method for improving model efficiency would be to omit the wind-wave element from the model system, as iterative wave solvers are computationally expensive. The Parametric Wave Model (PARAM) was formulated by Boyd and Weaver (2021) to overcome this burden. A parametric wave solution for wave growth with shallow water adjustments is used to calculate the wave field in restricted estuaries where the waves are locally generated wind waves and there is limited propagation of long waves or swell from the ocean. Outside the restricted estuarine portion of the domain, a third generation is used to calculate wave growth and propagation. When coupled into the ADCIRC hydrodynamic model to create ADCparam, run time is substantially faster than ADCIRC coupled with SWAN. Computational efficiency is greatly increased compared to third-generation wave and hydrodynamic coupled models, while solution integrity is maintained across the basin.
With the same goal of increased computational efficiency, the Multistage Mini-Ensemble Modeling System (MMEMS) is a system developed using one-way nesting to produce high-resolution predictions in estuarine domains [28]. Conventional modeling of detailed coastal zones requires blending the node density far offshore and limiting timesteps to that needed to resolve the most complicated elements. This results in excess calculation of data in regions of the grid that are largely ignored during analysis. MMEMS splits the grid into two stages of models: the first, simulating the large-scale coarse mesh, and, then, the second, using the solutions to force inlet boundary conditions of the estuarine fine mesh. Forecast runtimes are reduced by as much as 80% when compared to conventional single-domain models [28].

Parametric Wave Solver
To overcome the computational cost of solving for waves with discretized partialdifferential equations, the novel Parametric model removes the finite difference scheme altogether. As PARAM is designed for restricted estuarine simulations, the assumption can be made that the waves are locally generated and fetch limited. Wave parameters at every node can be solved using empirical equations for fetch limited waves. The assumption is weakest under large fetches and nominal wind regimes but has been shown to be valid under storm conditions. The system of parametric equations is recast to have wave parameters be functions of wind fetch, 10 m wind velocity, water depth, and bottom slope. There is no dependency on wave values adjacent in space or time. Since fetch values and bottom slopes are constant at every node, they can be pre-calculated for a given mesh, making water depth and wind velocity the only variables updated every PARAM timestep. Pre-processing for PARAM includes the calculation and tabulating of the constants: 360 degrees (in two-degree increments) of effective fetch length, bottom slopes, and average depths along fetch-length for every mesh node.
Systems of equations available in the Parametric ensemble are defined by Boyd and Weaver (2021) and include four sets of empirically derived significant wave height and period formulae from their contributors: Shore Protection Manual (SPM), Equations (1) and (5), [6,30,31]; Coastal Engineering Manual (CEM), Equations (2) and (6), [32]; Sverdrup-Munk-Bretschneider (SMB), Equations (3) and (7), [33]; and Texel-Marsen-Arsloe (TMA), Equations (4) and (8) [34]. TMA implements a wave spectra method, while the other three formulations use significant wave methods [35]. The SMB method is rooted in the SPM formulation, but recasts wind speed to represent wind speed at the water surface [33]. Significant wave period from the Malhotra and Fonseca study (2007) remains unmodified from the SPM method. The SMB equations have performed most similarly to SWAN in controlled test cases, followed closely by TMA and SPM [35]. PARAM can be adjusted to use the full four-member ensemble, or a custom ensemble of select formulae. In the test analyzed for this study, only SMB and SPM formulations are included. The CEM formulation, Equations (2) and (6), has degraded solutions in preliminary testing, due to its deep-water assumption [35], but remains an option in the model for further studies in alternative domains. The PARAM solution for significant wave height and period can be the average of some or all equations currently available.
T s (SMB) = T s (SPM) T s (TMA) = 1 Propagation in an estuarine domain can include shallow, intermediate, and deep-water waves, classified according to their water depth to wavelength ratio (d/L) [36]. Aside from the CEM formulation, which assumes only deep-water conditions, the ensemble members consider the full range of water depth classifications in their derivations. Depth variations upwind of a node also influence the solution and are accounted for with an inverse distance weighted average of depths along the relevant fetch length [35]. Water depths further away from a node contribute less to the value than depths nearby. Bottom slopes in each direction, used only if wave breaking is present, are determined from the difference between node depth, and depths at a 50-m step upwind.
Effective fetch is computed at each wet node by checking water depths at fixed step sizes in every direction; the same steps used in computing weighted depths and bottom slopes. Once a dry node is detected, the number of steps leading up to land and the step size, typically fixed at 50 m, are applied to the haversine formula for calculating fetch with respect to Earth's curvature [35]. Fetch lengths in individual directions are weighted according to all other fetch lengths within a ±45-degree arc to generate an effective fetch length [35]. Effective fetch improves wave calculations by considering the effects of fetch lengths that pass between, or next to, a land mass.
The effective fetch lengths, corresponding average bottom slopes, and the depths over fetch length, are calculated in 180 wind directions (360 degrees at 2-degree intervals) for every mesh node by pre-processing functions. Once tabulated by the pre-processing function into three files, they can be recycled into any PARAM model using the same gridded domain. Deviations in fetch and upwind depth with respect to water elevation within a hydrodynamic model are not considered substantial enough to render fixed values as obsolete. Even with large changes in water elevation, changes in depth weighting and fetch lengths do not cause anomalies in wave calculations [35].
By neglecting time rate of change to reduce computational demand, the ensemble equations 1 through 8 are diagnostic rather than prognostic. It is assumed that the time taken for the waves to approach fully formed states is negligible, given the fetch limitations and an extreme wind field. To mitigate error from this assumption, physical wave processes cannot be ignored. Wave height solutions for each parametric ensemble member are modified according to shoaling and breaking, as well as bottom friction effects. A bottom friction decay factor is calculated using a constant Manning's n coefficient specified by the user [31,37]. Once the wave parameters are calculated, they are used to determine wave-driven forces in a similar manner to ADCIRC + SWAN [38]. The PARAM model does not address wave refraction, diffraction, or reflection, and cannot resolve complex domains where these wave phenomena are dominant.
Validation for the wave output alone has been previously conducted for PARAM without ADCIRC coupling. From those simulations, the parametric ensemble's significant wave heights averaged a 12.9% mean absolute error using SWAN as the control, and a −5.5% difference when compared to isolated ADCP data [35].

Multistage Mini-Ensemble Modeling System (MMEMS)
The MMEMS is a one-way nesting approach for modeling estuarine circulation. A high-resolution estuarine mesh is nested within a coarse parent mesh. A broader, openocean domain provides elevation boundary condition forcing at the inlets of the embedded estuary, an estuarine high-resolution mesh. Using sparse nodal spacing and longer timesteps allows the coarse mesh model to run with minimal wall-clock time, even with SWAN coupling. Coarse and fine models are run in sequence, such that the high-resolution estuarine mesh is given the full non-periodic elevation boundary condition file from the parent model. The purpose of the ocean-scale coarse mesh is to feed the boundaries of the focal estuary and provide an accurate prediction of coastal circulation immediately outside the estuary.
The same two mesh and boundary information files were used without modification for the entirety of the study, focusing on the Indian River Lagoon, Figure 2. A highresolution mesh of the Indian River Lagoon, Figure 3a, stretching from North Mosquito Lagoon to Jupiter Inlet, was used for the embedded model. The mesh, adapted from a grid generated by [28], was refined for stability across a wide range of model types. The coarse parent mesh, Figure 3b was a modified version of the east-coast 1995 (ec95) mesh that includes the Gulf of Mexico, Caribbean Sea, and northwest Atlantic Ocean [39]. Spacing between nodes of the triangular elements ranged from 10 m in the high-resolution estuarine domain, Figure 3a, to over 400 m in the coarse ec95 mesh, Figure 3b, with increased nodal density around regions with complex bathymetry. Nodes with positive elevations were used in place of land boundaries to allow for wetting and drying.   Collections of two-stage models were wrapped into a complete system inspired by the ADCIRC Surge Guidance System [40,41], automating model set up and splitting MMEMS simulations into three steps: hindcast, nowcast, and ensemble forecast, Figure 5. The hindcast included a single ADCIRC run on the parent coarse mesh and served to spin up tidal and meteorological conditions from cold start prior to the nowcast. The generated hot start file was used to launch the parent grid coarse model simulation in the nowcast stage, which then fed boundary conditions to its nested high-resolution estuarine model. The two-stage nowcast applied waves to both the coarse and fine domains and ramped up all model conditions to produce hot start files for each of the forecast members. The nowcast ran from 1-3 days before the current cycle to the start of the meteorological forecast cycle, exceeding the required spin up time for the Indian River Lagoon mesh, Figure 3a) [28]. The coarse and fine mesh stages of forecast members started on this cycle and ran for the length of their respective meteorological inputs, or about 3 days. Waves might be omitted from coarse mesh stages to improve run time, but at the cost of model solution during storm events [28]. Left alone and with ample disk space, continuous power and connectivity, the MMEMS would continuously model forecast cycles as they were picked up from the NOMAD Server. Currently, seven meteorological members are available for MMEMS forecasts: NAM (North American Mesoscale, [42]), mean GEFS (Global Ensemble Forecast System, [43]), mean SREF (Short Range Ensemble Forecast, [44]), +1 standard deviations for GEFS and SREF, and −1 standard deviations for GEFS and SREF. The model spin-up uses NAM data records. While this study focused on the forecast ensemble, MMEMS is also capable of hindcasting, and event simulations. Input from a tropical cyclone or storm can be downloaded to simulate surge for the duration of the inclement weather.

Methodology
Structuring of the PARAM-MMEMS system began with the creation of coupled model ADCIRC + PARAM (ADCparam). Previously, PARAM was an independent serial application written in C++, which proved to be a challenge to integrate with ADCIRC's FORTRAN source code when compiled for parallel processing. The parametric wave program was redrafted in FORTRAN and modified to function as a subroutine within the ADCIRC source code. Other modifications were made to ADCIRC's main timestep routine, preparation processes, and output writers, resulting in a novel iteration of ADCIRC: ADCparam.
Regardless of the source code framework, the pieces of the coupled system function the same as their parent applications. ADCIRC's generalized global wave continuity equation is unaltered, and all constituents are customarily read in.
Development of ADCparam involved the integration of PARAM into the main AD-CIRC loop as a simple subroutine, Figure 6, that also required a series of cascading edits throughout the network of ADCIRC source code. Similar to ADCIRC + SWAN, the values passed from the wave model to the hydrodynamic model are radiation stress gradients, which must be normalized with water density before being added to ADCIRC wind stress gradients [45]. The PARAM subroutine was written into the primary timestep script to be prompted at a set wave calculation interval (RSTIMINC) if NWS = 5XX in the ADCIRC control file, where NWS was the specifier determining the form of wind and/or wave stresses to apply to the model. The interval RSTIMINC, typically set to 600 s, determined how often wave parameters and stresses were updated in the model. ADCIRC script reading in the control file was expanded to interpret the NWS specifier for PARAM waves, the wave calculation interval, and settings for time series output of wave heights, periods, and radiation stress gradients. Variables passed between the PARAM subroutine, radiation stress gradient function subroutine, main ADCIRC loop, and the output writer were all pre-allocated in a global script. Source codes initiating model conditions and output files were also amended to include PARAM parameters. In order to run ADCparam in parallel with multiple processors, prep functions were edited to partition PARAM input files for fetch, depth, and slope to match domain mesh file partitioning. All PARAM additions to ADCIRC source code included related parallel execution specifiers to note if indexes were in line with respect to local sub-processing loops. Within the PARAM function itself, the weighted averaging of depths along fetch lengths was altered in order to accommodate dry nodes upwind. If the averaged depth combined with nodal water elevation was negative, the true depth at node was used instead of the fetch length averaged depth. Similarly, the interpolation for wave heights and periods between adjacent 2-degree wind angles was also corrected. If either value implemented a negative water depth with its respective wind direction integer, the node was set to "dry" and wave values were null. Such instances occurred at nodes along the active water line, where wave-driven forces were imperceptible.
Editing the MMEMS to execute ADCparam required minor changes to configuration shell scripts, and estuarine fine mesh stage execution shell scripts, to link input files and control settings required for ADCparam. ADCparam input files and executables were added to MMEMS libraries. SWAN was still applied to the coarse parent grid model stage in MMEMS + ADCparam but was replaced by PARAM in the high-resolution estuarine stages. Aside from ADCparam additions, the only necessary update to MMEMS was on shell scripts downloading GEFS forecasts from NOAA databases. Following a recent change in NOAA server data restrictions, MMEMS had to download 3-h meteorology forecast files for the full cycle in sequence, instead of downloading all 24 files at once. Download rate restrictions only applied to the GEFS database, and did not affect download time of the other ensemble members.
A series of pilot analyses were conducted to refine the early versions of ADCparam. The initial programs ran on a closed version of the Indian River Lagoon fine mesh, Figure 3a, with only meteorological input. In one case, a constant 30 m/s wind coming directly from the North was applied, sustained for the entirety of a 3-day model duration and with zero ramp up. While this was unrealistic, it forced a maximum wind wave production under a severe-case storm surge scenario, particularly because a northerly wind allowed for peak fetch lengths across the domain. The pilot renditions of ADCparam were refined through these conditions, along with a simple MMEMS sub-model frame using NAM meteorological input, forced boundary conditions, varying friction, and geoid offset. Subjected to both wind and tidal frequencies, ADCparam evolved until operation was fully stable and output was satisfactory. An assortment of parametric ensemble combinations and Manning's n friction coefficients were tested throughout the iterative process, but no changes were made to hydrodynamic or core wave-resolving equations.

Model Setup
In the latest test to be analyzed, MMEMS + ADCparam forecasting was run for one full cycle, then duplicated at estuarine stages with original MMEMS + SWAN set up to compare output and runtimes. To evaluate MMEMS + ADCparam under the storm conditions that it is intended for, the real time test forecast aligned with the arrival of outer bands from tropical storm Elsa, Figure 7. Elsa was the first Atlantic hurricane of the 2021 season and decayed into a tropical storm before approaching the Straits of Florida [46]. As an early season storm, Elsa was also the first forecast event to apply a moderate wind regime to the Indian River Lagoon since MMEMS + ADCparam was first developed in Spring, 2021. While most of the Indian River Lagoon was not subjected to tropical storm winds (wind speeds over 17 m/s) and storm surge was minimal, the event was an opportunity to test MMEMS + ADCparam live forecasting, as well as to study the impact of PARAM replacing SWAN. The MMEMS + ADCparam test case used to exhibit real time forecasting was initiated at 1000 EDT 4 July 2021. This set the first forecast cycle period to start at 1200 EDT 4 July 2021, and end at 1200 EDT 7 July 2021. MMEMS orchestrated models with all 7 forecast ensemble members using meteorological and tidal input downloaded from public databases [28]. Since the forecast ran in real time using input from atmospheric forecasting products, MMEMS forecast members did not represent conditions from the actual track of tropical storm Elsa, which moved north of that forecast by the 4th July advisory, Figure 7, parallel to the West coast of Florida and made landfall near Apalachee Bay [46]. Meteorology forecasts sampled by MMEMS aligned with the predicted storm track released as part of the 4 July advisory, Figure 7.
To ensure replicability and draw direct comparisons between PARAM and SWAN, input files generated by the 3-day MMEMS + ADCparam real time forecast cycle were used to generate duplicate models of the estuarine forecast ensemble in both ADCIRC + SWAN and ADCparam. Replica models used identical input to the original real time forecast estuarine stages, but with a cold start spin up instead of building on hot start generation files from MMEMS nowcast members. Ramped up time in cold start simulations did not disrupt wave or hydrodynamic solutions within the time duration analyzed, as initial winds from outer storm bands first arrived 24 h into the forecast.
The MMEMS + ADCparam ensemble calculated forecasts with similar inputs used by Taeb and Weaver (2019). The system ran barotropic 2-Dimensional Depth-Integrated models in parallel with 20 allocated processors, with all 7 meteorological ensemble members (NAM, GEFSmean, SREFmean, GEFSmean + 1std, GEFSmean-1std, SREFmean + 1std, SREFmean-1std). Timesteps were 20 s for parent grid stages, and 2 s for the estuarine high-resolution stages. Since parent grid stages would be identical, only the estuarine high-resolution stages were duplicated with ADCIRC + SWAN and ADCparam cold starts, using the same 2 s time steps and 20 processor partitions. In the SWAN and PARAM scenarios, waves were calculated every 10 min. All input was identical, aside from files needed for the respective wave solvers, including the tidal inflow at boundary nodes from the parent grid, output specifications, and nodal attributes [28]. Meteorological input in all ensemble members were read by ADCparam and ADCIRC + SWAN in the same format with 3-h timesteps. A nodal attributes file specified spatially varying friction and roughness coefficients, sea surface height above geoid to account for seasonal water levels, surface canopy and roughness coefficients to reduce or turn off wind stresses over land, and primitive weightings in the generalized wave continuity equation.
For this case study, the PARAM wave solver was set to use the SMB and SPM ensemble members, Equations (1), (3), (5), and (7), based on model performance in preliminary testing and in analysis by Boyd and Weaver (2021). SMB and SPM ensemble members used the same methodologies to calculate wave height and utilized the same parametric equation for significant wave period, Equation (5).
The PARAM subroutine used a constant Manning's n friction coefficient, rather than assigning values unique to each node, as specified by the ADCIRC nodal attributes file. A separate preprocessing code was written to read ADCIRC grid and nodal attributes files, and, then, a single depth-weighted average Manning's n was computed to represent the full IRL wetted domain. The calculated Manning's n coefficient of 0.027 was used in PARAM formulae for physical processes, modifying the significant wave height of each Parametric wave ensemble member. In contrast, the SWAN wave solver was passed spatially varying Manning's n friction coefficients from the ADCIRC nodal attributes file. Future iterations of the PARAM model will include this feature.

In-Situ Measurements
Water elevations averaged across all sets of estuarine stage forecast models were evaluated with data from in-situ measurements from Onset Hobo U20-001-01 pressure sensors [47]. Hobo sensors were deployed at the indicated gauge locations, Figure 8, for cycles lasting 1-3 months. The instruments were deployed by mounting them at the base of pilings under approved fishing and boat piers. While deployed, the sensors recorded temperature, water pressure, barometric pressure, and computed water levels in 5-min intervals.

Results
Run times were recorded for every individual model, and global output files for wave heights, water elevation, and wave periods were visualized in SMS 11.1 (Aquevo, Provo, UT, USA) [48]. Timeseries output was analyzed at select stations chosen for their varying sensitivities to tidal flow and fetch in North/South directions, Figure 8.
Statistical analysis between ADCIRC + SWAN and ADCparam model output used average bias, Equations (10) and (11) and Normalized Mean Absolute Error (NMAE), Equation (12), averaged across all nodes at one timestep, all timesteps in the sample wind event duration (from the 1.25-day mark through the end of the 3-day cycle) at a single node, and all nodes across all timesteps in the sample duration. A similar normalization was applied to model run time reduction, Equation (13). Average bias was also applied to compare both models against in-situ pressure gauge data, Equation (11).

Model Runtime
The main goal of running the test was to compare run times between ADCparam and ADCIRC + SWAN, demonstrating the purpose in developing the coupled Parametric model. Through this study we determined that the newly developed ADCparam reduced run time from ADCIRC + SWAN for individual models by an average of 51%, Equation (13), Table 1. Run time reduction was compounded with every fine grid model run in the MMEMS + ADCparam system, such that when all 7 members of the ensemble were run within the real time forecast 3-day cycle, 733 min of wall time were saved, Table 2. This was an overall reduction in run time from MMEMS with ADCIRC + SWAN of 33%, Equation (13), after also including run times from parent grid stages of MMEMS. Total MMEMS model run times for a single 3-day cycle combined wall times of one hindcast coarse parent grid model, one coarse parent grid nowcast model, one estuarine fine grid nowcast model, seven coarse parent grid forecast models, and seven estuarine fine grid forecast models. Both MMEMS with ADCIRC + SWAN and MMEMS + ADCparam used identical parent grid models with SWAN in the nowcast and forecast stages, and no waves in the hindcast stage, Figure 5. Hindcast run times were negligible, but parent grid models of nowcast and forecast stages lasted 80 to 90 min each when applying SWAN (running in parallel with 20 processors).

Model Output
Model impact was determined by comparing the water elevations, significant wave heights, and spectral peak wave periods between ADCparam and ADCIRC + SWAN. Figure 9 provides insight into the meteorological forcing of each ensemble member, and the influence on wave input parameters (effective fetches and inverse-distance weighted average depths). From day 1.25 through to the end of the model, wind velocities increased and had sustained direction, indicating the arrival of outer bands from Elsa, Figure 9a,b. Wind speed magnitudes were nominal in the first 24 h, and SWAN solutions ramped up from zero in the first 6 h. ADCIRC + SWAN required no more than a 24-h ramping duration for model spin up, [28]. The time frame of higher wind speeds, the 60-h timestep in the 3-day cycle, representing conditions on the periphery of tropical storm Elsa, was analyzed for bias and NMAE values averaged over time, Figure 9b. Omitting the initial 1.25 days, the mean bias values for the model (averaged across all nodes at all subsequent timesteps) were: −0.189 mm water elevation, −0.0303 m wave height, and 0.0631 s wave period. Corresponding average NMAE values were 0.23% for water level, 51.53% for wave height, and 10.44% for wave period, Table 3.  Figure 9a,b. After the first day, wind direction within individual ensemble members was generally consistent, Figure 9a. Low variation in wind direction reduced duration limited wave growth, increasing the accuracy of fetch limited solutions from ADCparam. The dominant wind direction forecasted by the ensemble was out of the Southeast, reinforcing contour gradients of both model solutions in Figure 10. Contours were centered at −0.36 m due to model geoid offset, Figure 10a,b. For the chosen input conditions, water level was nearly identical. At the 60-h timestep of the model, Figure 10, the ADCparam, Figure 10a, and ADCIRC + SWAN, Figure 10b, solutions differed by less than 1%, Figure 10c. The average water level bias, Equation (10), across all nodes was −0.26 mm. Given the moderate wind field, minimal variations were expected.
In areas nearest the inlet, where tidal influence was greatest, the elevation differences between the individual ensemble solutions were negligible, Figure 11a. The Sebastian Inlet station had the greatest error of the four nodes, with an average bias of 1.7 mm and NMAE of 0.42%. Larger deviations between model systems were most pronounced at locations where the tidal amplitude was small, as in Figure 11c. Set up and set down from winds also became more noticeable with less tidal influence, Figure 11b-d. At all three locations, the range of forecast model results reflected similar phases and amplitudes to in-situ pressure data. Near Lansing Island and Kars Park, where tidal influence was minimal, pressure readings largely remained within the bounds of the minimum and maximum water levels predicted by the full ensembles, Figure 12b,c. Since total variation in water level was small near Lansing Island and Kars Park, common trends in set up, Figure 12c, and set down, Figure 12b, were visible. Deviation was largest in Sebastian, where the station was subject to larger tidal amplitudes, Figure 12a. Noise in the pressure gauge signal was present at all three stations. Despite applying a 10-point (50 min) moving average to pressure gauge data sampled every 5 min, noise was still noticeable at Kars Park where overall changes in water level were smallest. ADCIRC + SWAN and ADCparam demonstrated near-identical error compared to station data. Average bias from in-situ data over the time series was 3.2 cm at Lansing Island, Figure 12b, for both ADCIRC + SWAN and ADCparam models, Equation (11). At Sebastian, Figure 12a, average bias was 4.8 cm and 4.7 cm for ADCparam and ADCIRC + SWAN, respectively. Bias values were near −0.8 cm for ADCparam and −0.7 cm for ADCIRC + SWAN at Kars Park, Figure 12c.
While significant wave height would be the preferred measurement to verify the ADCparam model, readily available water elevations recordings began to show the impact of replacing SWAN with PARAM. The PARAM hydrodynamic contribution induced marginally higher magnitudes for set up and set down than SWAN because of the more conservative wave height solution provided by PARAM. The impact of replacing SWAN with PARAM on storm surge predictions would be subtle, as both sets of ensemble model were equally incorrect. Error in both models was primarily attributed to local wind variations that were not captured by the low spatial and temporal resolution of all 7 meteorological ensemble inputs, as well as discrepancies in tidal forcing against historical tidal signals [49]. The most spatially refined atmospheric model used, NAM, had a resolution of 12 km, and the shortest temporal resolutions, applied in SREF ensembles, were 3-h cycles [42,44].
Modeled significant wave heights for ADCparam and ADCIRC + SWAN revealed similar patterns of wave height distributions. Compared to the ADCIRC + SWAN results, the ADCparam predicted a larger significant wave height in shallow and restricted regions and slightly smaller wave heights in the most open regions of the domain.
Due to the fact that the SWAN solution is time dependent, there was a contrasting spin up and slight phase lag in the SWAN solutions, Figure 13, compared to the PARAM subroutine. These differences were expected, and the overall offset in amplitude was marginal. The ADCIRC + SWAN model was given ample time to spin up before winds increased in magnitude at the 24-h mark. In the Northern area of the domain, Figure 13b,c, the ADCIRC + SWAN solution generally predicted more extreme maximum and minimum wave heights than ADCparam, although both produced a broad ensemble spread. The mean solutions were most in agreement north of NASA causeway, Figure 13c, with only 4.75% NMAE. In most cases, the ADCparam ensemble was overpredicting wave heights compared to SWAN, as expected, Figure 13a,d. In focal areas where waves developed over a long fetch, SWAN wave heights began to eclipse PARAM wave heights. Domain-wide wave height bias and NMAE at the 60-h timestep were −0.029 m and 43.75%, respectively. The wave heights near Sebastian Inlet showed the largest difference between the two models, Figure 13a, with an NMAE of 53.80% and an average bias of −0.044 m (ensemble mean at individual station node averaged over all timesteps after the 1.25-day mark). Effective fetch lengths and average depths were smaller at this location, Figure 9(c1,d1), generating smaller waves that contributed less to storm surge.
Peak wave period solutions featured a less pronounced contrast in magnitude between SWAN and PARAM than wave height output, and the trend in periods was similar between the two models. SWAN predicted a slightly larger spectral peak period, but remained nearly in phase with PARAM, Figure 14. Compared to ADCIRC + SWAN, ADCparam slightly underpredicted the wave period predicted by SWAN across most of the domain. Average bias and NMAE at the 60-h timestep were 0.09 s and 11.34%, respectively. The largest discrepancies in Figure 14 were 10.37% NMAE and −0.12 s average bias near Sebastian Inlet, Figure 14a, (ensemble mean at individual station node averaged over all timesteps after the 1.25-day mark). Similar to the wave height solution, ADCparam predicted a larger wave period than ADCIRC + SWAN in nearshore regions with smaller depths and fetch lengths, where waves were beginning to develop, Figure 14a. Across the rest of the domain, where waves were more mature, the models aligned more closely. Near Pineda Causeway, mean NMAE was 2.45%. Here the profiles appeared most in phase and shared a similar maximum and minimum, Figure 14b. Minor phase lag across all stations stemmed from PARAM omitting time dependent wave growth. Melbourne causeway and NASA causeway nodes yielded the most pronounced discrepancies in ensemble range between maximum and minimum, Figure 14c,d. The SWAN solution for period displayed subtle signs of instability in oscillating sections of the time series near the Melbourne and Pineda causeways, Figure 14b,d, as well as where the smallest ensemble member returns zero wave periods at the 1.5-day mark near NASA Causeway, Figure 14c. Unstable oscillations in SWAN wave period near Pineda and Melbourne Causeways, Figure 14b,d, were caused by abrupt reductions in GEFSmean and GEFSmean-1std wind speeds, the minimum ensemble members in Figure 9(a2,a4). When wind speeds approached zero in the GEFSmean member, the ensemble minimum in Figure 9(a3), PARAM wave heights and periods approached zero and SWAN wave heights and periods were at zero, Figures 13c and 14c.
The rise in wind speeds across all ensemble members coincided with increases in wave height and period, Figures 13 and 14. As the ensemble had a narrower field of wind speeds near Sebastian Inlet, Figure 9(b1), the ranges of model solutions for wave heights and periods were more dense, Figures 13a and 14a. Ensemble forecasts encompassed a diverse range of effective fetch lengths, given wind directions of respective members, but the longest fetches after the 1-day mark were near Melbourne Causeway, Figure 9(c4). Waves from south-southeast winds were allowed to propagate unobstructed for over 12 km in this region of the Indian River Lagoon, resulting in larger wave heights and periods, Figures 13 and 14. The Sebastian Inlet and Pineda Causeway areas, Figure 9(c1,c2), had shorter fetches, and, consequently, smaller wave heights and periods.
Ensemble ranges for inverse weighted average depths at each location, Figure 9d, followed the same trends as fetch length, Figure 9c. As winds travelled over a longer stretch of the Indian River Lagoon, they transverse over deeper water, increasing the average depth along the line of fetch, Figure 9(d3,d4). Regions leeward of a barrier island were represented with shallower depths averaged over shorter fetches, Figure 9(d1,d2). As with long fetch, larger average depths induced greater PARAM wave heights and periods that aligned more closely with SWAN solutions, Figures 13 and 14.
While high frequencies were present in the ramping durations and regions of most wave models, wave periods lower than 0.7 s were outside the scope of SWAN design solution settings. Peak periods and significant wave heights were integrated from the computed spectrum over a frequency range with a user defined minimum, 0.03 Hz, and maximum, 1.43 Hz. Within the weaker ensemble members, durations of light wind forcing did not produce waves within the design range and compromise sections of the SWAN solution, as indicated by instabilities along the lower boundary of ADCIRC + SWAN wave period ensemble spread, Figure 14. This contributed to regions of larger NMAE in restricted regions, Figure 14c. Across most of the domain, the ensemble average wave solutions were largely within the constraints of integration, and agreement between PARAM and SWAN suggested uncompromised SWAN computations, Figure 14.
Average bias values in wave period for the full ensemble, Table 3, as well as for selected stations and contours, Figure 14, were slightly positive. SWAN calculated a larger wave period than PARAM at nodes with the largest fetch lengths and most developed waves, implying a minor difference in wave spectra between the two models. SWAN solutions were based on a default JONSWAP spectrum, which has been verified extensively in open ocean environments, such as the North Sea and northeast Pacific Ocean, and is not a direct function of depth, [50,51]. The PARAM empirical formula for spectral peak wave period, Equation (5), was a depth-dependent equation that was developed for shallow water by accommodating energy loss due to friction and percolation, [6,52].
In areas with large average depths and effective fetch lengths, PARAM error in wave height was low compared to SWAN, Figure 13c,d. SWAN wave periods grew slightly larger than PARAM as depth and fetch increased, but NMAE remained lower than NMAE, when comparing significant wave height, Figure 14c,d.
Individual ensemble members showed the same trends between solutions and windwave parameters, Figures A2-A4. PARAM wave heights were in closer agreement compared to SWAN in ensemble members with larger sustained fetch and depth, such as GEFSmean-1std near Melbourne Causeway, with an NMAE of only 10.09%, Figure A3d. At the same node near Melbourne Causeway, GEFSmean + 1std had the most error with 81.86% NMAE, which was most likely attributed to the nominal effective fetch, Figure 9(c4), and lower depth, Figure 9(d4), given the wind direction of that meteorological forcing. ADCparam ensemble members had lower differences in wave period when compared to ADCIRC + SWAN in areas where effective fetches were moderate. Wave period NMAE at the stations ranged from 3% at Pineda causeway in the SREFmean ensemble member, to 24% at Sebastian Inlet in the GEFSmean ensemble member, Figure A4a,b. All four wave period time series output stations had NMAE under 8% within the NAM ensemble member, Figure A4. Average wave period bias was often negative where fetches and depths were small, and slightly positive where fetches and depths were larger.

Discussion
Qualitatively, the model output shows that ADCparam provides a hydrodynamic solution nearly equal to ADCIRC + SWAN, wave parameters are well within the range of expectations, and the wave height and period solutions are on the order of what is predicted by SWAN. Water elevation output between ADCparam and ADCIRC + SWAN is similar, with discontinuities mostly around complex boundaries and grid segments. There are no notable errors in water level, and comparison to in-situ measurements supports the model's credibility in both its forecasted tidal signals, and its wind-driven hydrodynamic processes. Discrepancies between the numerical models and in-situ measurements were attributed to ADCIRC + SWAN model sensitivity to spatial resolution of meteorological forcing [49], as wind and pressure fields only reached 16 km 2 resolution for the most refined ensemble. Noise was also prominent in the measured data due to the HOBO sensor sampling system. Each water level sample was recorded at an instant in time instead of across a set of bursts, allowing free surface displacement from waves to distort the water elevation signal. No large discrepancies between velocity field outputs from the two models were observed during the study, and forecasted current velocities could be verified with data from ADCP instruments deployed in the Indian River Lagoon during future storm events.
The average NMAE and bias across the full model agreed with values calculated at select nodes and timesteps. Each of the time series figures highlighted the contrast during SWANs initial ramping. Overall, the error in water elevation was close to zero. SWAN wave periods were on par with PARAM for most of the model, and NMAE confirmed that they were proportionately similar. Wave heights exhibited the same correlation, but with larger NMAE and a negative bias across most of the domain. Station output, Figure 13, and corresponding wind parameters, Figure 9, supported expectations that wave heights between the two models were more alike in deeper, open-fetch regions of the basin. Discrepancies where wave heights were small and less developed, such as the leeward side of Merritt Island, amplified the overall NMAE in wave height.
The discrepancies in wave height and period showed the contrast between parametric calculations and the SWAN implicit solver. Wave parameters give a detailed picture of how each method behaves differently and indicate that SWAN and PARAM would be in agreement over most of the domain if wave periods were not slightly underpredicted in PARAM (or overpredicted by SWAN), and wave heights were not overpredicted by PARAM (or underpredicted by SWAN). In domain regions where wave height and period magnitudes for SWAN and PARAM were offset to various degrees, shape of time series plot profiles were still in phase. Variation in wave height and period solutions of smaller magnitudes, as well as minor oscillations in SWAN, Figure 14d, might also be a result of durations of nominal wind regime applied to the individual models within the ensemble, and presence of changing wind directions. Frequency limiters within SWAN and challenges in resolving short-period waves further emphasized the value in an unconditionally stable parametric wave model, such as ADCparam.
SWAN and PARAM showed better agreement under stronger wind forcing, and moderately open regions of fetch and depth. The parametric ensemble approached an upper fetch limit for accurate wave period solutions, where SWAN continued to generate longer waves as fetch increased, Figure 14. PARAM wave height computations were more divergent than wave period in shallow, restrained regions, but had a higher threshold of fetch length before agreement with SWAN began to decay. The contours, along with relations between wind parameters, Figure 9, and wave solution time series, suggested a fetch "sweet spot" with ideal fetch, depth, and wind magnitude for the parametric ensemble equations used. While performance might change under critical wind conditions, specific members of the parametric wave solving ensemble potentially skewed solutions and could be modified or removed to better fit the wave magnitudes of SWAN over a wider range of averaged depths and effective fetch lengths. The results shown exhibit an offset in wave spectra, where SWAN's default JONSWAP spectrum and upper frequency limits shift solutions towards longer period waves. More in-situ data is required to determine which model best represents the wave climate in an enclosed shallow estuary, and which parametric ensemble combination within PARAM is most accurate.
The function of ADCparam was successfully demonstrated by the substantial increase in computational efficiency. SWAN run time is entirely dependent on its number of iterations, causing wall-clock time differences to fluctuate, based on model parameters. Further tests with varied meteorological input and mesh domains will pinpoint the true range of reduced wall-time as compared to the most efficient model setup for ADCIRC + SWAN. SWAN might also be omitted from nowcast and forecast coarse mesh model stages, Figure 5, further reducing run time for both MMEMS + ADCparam and the original MMEMS. Walltime reduction compared to single domain ADCIRC models was estimated at 89% based on prior model analysis, as the MMEMS already reduced run time from single domain models by as much as 80% [28]. From this ensemble study, ADCparam wall clock time was a moderate improvement, 51%, over ADCIRC + SWAN, leading to overall runtime improvement of 33% over MEMMS + SWAN.

Conclusions
Coupling the PARAM into the existing ADCIRC framework to create a parametric wave solver option for restricted estuaries can improve the efficiency of coupled model simulations by reducing the wall-clock time by over 50 percent. This pilot study confirms the potential of PARAM as a valuable addition to MMEMS and ocean circulation modeling. ADCparam run time is 51% faster than an ADCIRC + SWAN model, and hydrodynamic solutions align, given the moderate meteorological regime.
ADCIRC + SWAN was selected as the traditional model for evaluating ADCparam based on aptitude of SWAN for application in shallow coastal environments, and past MMEMS systems using ADCIRC + SWAN in estuarine stage models [28]. SWAN is representative of third-generation wave models in this study, but further analysis of PARAM will include a wider base of phase averaging wave models. The core function of PARAM is also applicable to most hydrodynamic models, and should be tested within packages other than ADCIRC, such as Delft3D or SLOSH [12,13]. The parametric subroutine is structured such that it can be expanded to include additional ensemble members. Currently, the wave ensemble can use adapted versions of SPM, SMB, CEM, and TMA equations. When each of these were individually compared to SWAN in previous studies, the SMB and TMA schemes outperformed CEM and SPM, with CEM appearing to degrade the ensemble average [35]. However, in less critical wind regimes, such as the meteorology applied to this study, the smaller wave heights generated by SPM and SMB formulae (selected for the Elsa case study) mitigate bias associated with assuming fully formed fetch limited waves. Forcing a major storm surge event would further align with the PARAM assumption of fully formed fetch limited waves and justify the use of more conservative ensemble members. This conclusion is dependent on physical processes included in the model and implies that the wave equations are suited to different conditions according to the empirical data they are derived from. Select wave ensemble members will perform better under alternative domains and storm conditions. For instances where the performance of each member is dependent on the water depth classification (depth/wavelength), then weighted averaging can be applied accordingly. Additions have been considered for the expansion of the Parametric ensemble, including the Donelan method and the Krylov method. The Donelan method, in particular, is ideal for fetch limited environments, as it was developed for the Great Lakes [23].
Wave solutions will improve as the Parametric model is further developed and tested, but even after the PARAM model has been optimized, it is not meant to fully replace third-generation wave models. It will fulfill its role as an economic means to include waves in estuarine hydrodynamic model simulations, while iterative wave solvers remain crucial for large-scale ocean models tracking the effects of more significant swell. The larger course mesh of MMEMS will continue to use SWAN accordingly. Within PARAM's niche of fetch limited environments, inclusion of waves in a numerical model no longer needs to be a computationally intensive undertaking.   Figure A1. Each of these locations is selected based on location in the domain and associated fetches (station locations are specified in Figure 8).  Figure A1. Each of these locations is selected based on location in the domain and associated fetches (station locations are specified in Figure 8).  Figure A1. Each of these locations is selected based on location in the domain and associated fetches (station locations are specified in Figure 8).