Development of a Flexible Coupling Framework for Coastal Inundation Studies
Saeed Moghimi Andre van der Westhuysen Ali Abdolali Edward Myers Sergey Vinogradov Zaizhong Ma Fei Liu Avichal Mehra Nicole Kurkowski
physics.ao-ph















address1,address2]Saeed Moghimi mycorrespondingauthor
[mycorrespondingauthor]Corresponding author
(Saeed Moghimi; Email: Saeed.Moghimi@noaa.gov)



address3]Andre Van der Westhuysen _andre
address3]Ali Abdolali _andre
address1]Edward Myers _andre
address1]Sergey Vinogradov _andre
address3]Zaizhong Ma  _andre
address4]Fei Liu  _andre
address3]Avichal Mehra _andre
address5]Nicole Kurkowski _andre

[address1]​NOAA Coast Survey Development Laboratory, National Ocean Service, Silver Spring, USA. 
[address2]University Corporation for Atmospheric Research, Boulder, USA. 
[address3]​NOAA Center for Weather and Climate Prediction, College Park, USA. 
[address4]​NOAA Earth System Research Laboratory, Boulder, USA. 
[address5]​NOAA Office of Science and Technology Integration, Silver Spring, USA. 

[_andre]Co-authors contribution:
 
1. Saeed Moghimi: Responsible for development of the ADCIRC, WW3-data
and HWRF-data Caps and developing the coupled application. Writing the body of the manuscript and plotting figures.

2. Andre Van der Westhuysen: Responsible for unstructured WW3 cap
development.
Inclusion of the WW3 cap in the coupled application and running full coupled cases; manuscript preparation and analysis of the results.

3. Ali Abdolali: Responsible for WW3 cap development and running
stand alone wave model, plotting wind and wave related parts of manuscripts, manuscript preparation and analysis of the results.

4. Edward Myers: Scientific support, discussion and project
management; 
5. Sergey Vinogradov: Scientific support and discussions.

6. Zaizhong Ma: Responsible for wind field modeling using HWRF.

7. Fei Liu: Support for ESMF/NUOPC cap development.

8. Avichal Mehra: Scientific support and project management.

9. Nicole Kurkowski: Scientific support and project management.




To enable flexible model coupling in coastal inundation studies, a coupling
framework based on ESMF/NUOPC technology under a common modeling framework
called the NOAA Environmental Modeling System (NEMS) was developed. The
framework is essentially a `software wrapper around atmospheric, wave and storm 
surge models that enables its components communicate seamlessly, and efficiently 
run in massively parallel environments.   
We  implemented  the coupled application  including  ADCIRC and
unstructured WWAVEWATCHIII caps as well as NUOPC compliant caps to read Hurricane Weather Research and 
Forecasting Model (HWRF) generated forcing fields.    
We validated the coupled application for a laboratory test and a full scale
inundation case of the Hurricane Ike, 2008, on a high resolution mesh covering the whole US
Atlantic coast.
We showed that how nonlinear interaction between surface waves and total water level
 results in significant enhancements and progression of the inundation
 and wave action into land in and around the hurricane landfall region. We also
 presented that how the maximum wave setup and maximum surge regions may happen
 at the various time and locations depending on the storm track and geographical
 properties of the landfall area.




wave-current interactiontidal inlet river plume wind waves three-dimensional circulation  wave breaking 









§ INTRODUCTION








To establish a coastal flooding modeling system, several model
components based on the target geographical region need to be coupled.  To
accurately simulate the total water level in a tropical hurricane land-falling
inundation study, a dynamically coupled system of numerical models including
storm surge, surface waves, inland river flooding and numerical weather
prediction are necessary. On top of that based on the geographical location
other model components may need to be employed. For instance, to setup an
efficient coastal flooding prediction system for Alaska region, inclusion of a
sea-ice model is essential. 






In recent years, Earth System Models were proven to be invaluable tools that
enabled us to better understand and more accurately predict our environment. 
Each system includes a coupled applications that consists of several model
components to represent relevant physical processes.  The model components
are expected to interact with each other similar to what takes place in nature.

There are several Earth System Model software flavors that enable model
components to communicate by importing and exporting data
<cit.>.
The Earth System Modelling Framework (ESMF) has been utilized to develop
several earth system coupled applications worldwide
<cit.>. To
increase ESMF interoperability, the National Unified Operational Prediction Capability 
(NUOPC) consortium developed a layer consisting of a set of generic components
 <cit.>. 
NUOPC layer is a software wrap around ESMF and was developed
collaboratively by several research and operational centers. The primary
objectives behind NUOPC design are to be reusable, extensible and
portable framework for ESM coupling.





The NOAA Environmental Modeling System (NEMS) is a coupled modeling
infrastructure designed to address increasing needs for prediction of the earth
environment at a range of time scales.  NEMS includes several external
model components that have a primary source code outside NOAA. Therefore, NOAA
only needs to maintain and develop the coupling interfaces (so-called model
caps) of the given modeling component. In turn, the NEMS ecosystem allows
connecting various combinations of model components into a number of different
coupled model applications to address specific environmental phenomena at
specific time scales.         

The present research goals are to develop a flexible and generic coupling
between ADvanced CIRCulation model <cit.>
and WAVEWATCH III <cit.> via their respective
NUOPC caps, and to provide an infrastructure to make future development and inclusion of various model
components, such as river and inland flooding coupling, seamlessly possible. The
current development of the NUOPC caps provides the possibility to perform
dynamical coupling of ADCIRC and the unstructured version of WW3, as
well as various atmospheric models ATM. The cap developed for ADCIRC
is capable of importing atmospheric forcing and surface wave fields, and exporting water
surface elevation and current velocity to the connected model components.
Conversely, the cap developed for unstructured WW3 imports atmospheric
forcing, water levels and current, and exports the wave radiation stresses
required to force ADCIRC.         

The first application of this new coupled system is the so-called Named Storm
Event Model (NSEM), a high-fidelity model for hindcasting coastal
inundation and total water level. It is being developed to meet the requirements
of the NOAA's Consumer Option for an Alternative System to Allocate Losses
(COASTAL) Act of 2012. This modeling system includes ADCIRC as the 
hydrodynamic component, WW3 as the wave model, the Hurricane
Weather Research and Forecasting Model (HWRF) as the atmospheric component
<cit.>, and in future the National Water Model
(NWM) as the inland  hydrological component <cit.>, see Fig.
<ref>. 


The structure of this paper is as follows. First, we describe the envisioned
design of the NEMS ADCIRC-WW3 coupled application and the methodology. 
Then a detailed description of the ADCIRC cap implementation and available
coupling options is given. This is followed by a similar description of the
WW3 cap. Subsequently, we present the results of the coupled system. Finally 
we present verification of the coupled ADCIRC-WW3 application for the
laboratory flume case of <cit.>, as well as a full scale
storm surge inundation event during Hurricane Ike, 2008 in the US Gulf of
Mexico.



§ STRUCTURE OF THE COUPLED APPLICATION



A typical NUOPC application includes a number of generic components that
provide an interface to the underlying ESMF infrastructure for generating
and operating a coupled application in a fairly straightforward and seamless
manner. The generic components are defined as follows:    

  
  * A Driver manages all the components to initialize, run,
  finalize and keep track of time for exchanging information among model
  components.    
  
  * Connectors are used to execute field matching, grid remapping
  and data redistribution among model components.  
  
  * A Model (cap) wraps each model component code (e.g. ADCIRC
  and WW3) to provide  a generic interface and standard metadata suitable
  to be plugged into the Driver, and form a multi-model coupled
  application. 
  
  * An optional Mediator wraps custom coupling code to calculate
  quantities which  includes data from several model components or requires
  operations such as time averaging. 
 

The system includes methods and utilities for time management, error handling,
high performance inputs/outputs (I/O), grid remapping and field
interpolation.  Since NUOPC is a layer around ESMF library, function
calls to both NUOPC and ESMF are possible and sometimes are
necessary.       

In this research, we developed a NUOPC application that includes a driver,
three NUOPC enabled model components and four connectors. The components
are not allowed to directly access each other's data. The only way the data
moves in or out of a component is via instances of an ESMF_state class.
The state is a container that wraps native data and also includes a metadata to
let the other components know about name, coordinates and decomposition of the
actual packed data. 






The driver component accesses ADCIRC, WW3 and ATMesh model
components via their SetServices() methods. It reads basic information for
how to initialize and run the model components from a configuration text file 
(Fig. <ref>).  The configuration file contains information about
name of the model components, number of processes to be associated to each model
component, the coupling time intervals, and the order of data exchange among the
components.  The driver also initializes the number of connectors by providing
the name of the sending and receiving model components.  Therefore, for a
dynamical two-way coupling between two model components, two connectors are
required.           

The connector component initializes at the run time by matching the list of
available import and export fields advertised by the model components. The
connector establishes the connection based on matched import and export fields.
The connector also has access to the domain decomposition and computational
domain discretization of the connected model components. It will generate a
remapping and necessary weight matrices for interpolation of the fields among
model components at the initialization phase. In other words, the connector
receives exported data in the form of an ESMF_state in the native grid or
mesh from the exported model component and passes it to the importing model
components in its own native grid or mesh definition. This remapping facility
allows the coupled system the freedom to use different meshes for, say, the
circulation and wave modeling components, and/or different domain
decompositions. This is useful in cases where the wave model component requires
a different mesh optimization to resolve its distinct physics, or more
computational cores for load balancing.               

The ATMesh cap was developed as a placeholder interface for a full live
atmospheric model, which was not included in our NSEM application due to scope
limitations. This so-called data cap reads weather prediction outputs (from a
NetCDF data file), initialize required NUOPC/ESMF objects and provide
requested data and information to ADCIRC and WW3 caps via the
NOUPC/ESMF backbone. The NWM hydrological model
component and its associated connectors are not yet implemented (Fig. <ref>).       

 




§ ADCIRC MODEL


The ADvanced CIRCulation model (ADCIRC) is a finite element hydrodynamic
community model originally developed by <cit.>. ADCIRC 
is undergoing continuous development by groups of scientists and engineers. Its
natural finite element unstructured mesh capability, and several modules
specifically addressing various aspects of the coastal flooding and tropical
cyclone forcing, make it one of the best tools available for coastal inundation
studies. ADCIRC operates in either two-dimensional depth-integrated (2D)
depth-averaged (barotropic) and three-dimensional (baroclinic) modes.  In the 2D mode, it
solves equations for both water surface elevation and the depth-averaged
velocity fields.  For more details about ADCIRC governing equations,
numerical methods and wave forcing implementation please see
<cit.>.            

ADCIRC is written in modular FORTRAN and supports parallel execution
on massive supercomputers using MPI architecture.  The code structure is
partitioned in three distinct initializing, running and finalizing phases ready
for the ESMF coupling.  The model initializes by a call to ADCIRC_Init()
which also receives a MPI  communicator from the driver. The subroutine
reads necessary input files for constructing the computational mesh including
nodes location and connectivity.  It also builds a local and global nodal map to
reference which nodes reside on which MPI process, and to identify their
global relationships.  It reads input information to constrain the model such as
bathymetry, meteorological forcing, and freshwater inflow and open boundary
conditions.  As a part of the initialization, ADCIRC also checks and connects to
all requested output files that will be used as containers to fill in the model
results.             

ADCIRC enters the run phase by a call to ADCIRC_Run() subroutine,
which also receives an argument for the number of time steps (NTIME_STP)
for that specific run request. The start time and end time of the simulation is
determined during the initialization phase. The model run takes place via a time
loop in which, at every time step, a single call to the TIMESTEP()
subroutine occurs. All the computational steps for applying forcing and boundary
conditions to produce the final results are being performed in this subroutine.         
The ADCIRC concludes its run by a call to ADCIRC_Final() subroutine
where some of the final post-processing and check for MPI finalizing are
performed.    

 



 §.§ ADCIRC coupling interface (cap)




The ADCIRC NUOPC cap performs the coupling in all the three phases:
initialize, run and finalize. In the development of the NUOPC cap for
ADCIRC, extreme care and attention were paid to minimize changes to the
original ADCIRC code. At the initialization of the NUOPC application,
a global MPI communicator is created by ESMF infrastructure and a
dedicated set of processes passed to ADCIRC via a MPI communicator
based on the number of processes requested for ADCIRC in the configuration
file. At the initialization, ADCIRC cap also gets connected to available
import and export field matches accepted by the communicators.          


After information exchange among the model components, the
ModelAdvance() subroutine of the ADCIRC cap calls the
ADCIRC_Run() subroutine to perform the next run interval. Tab.
<ref> shows the list of the exported and imported fields
currently accepted by the ADCIRC cap.  The naming conventions of these
variables are defined in the NUOPC field dictionary to allow
interoperability with other NUOPC components. We modified and tested ADCIRC
preprocessing and main model code to accommodate various coupling arrangements. 
The NWS input parameters in fort.15 input file
are described in Tab.  <ref> <cit.>.          







§ WAVEWATCH III


WAVEWATCH III (WW3) <cit.> is a third-generation spectral wave model
that solves the wave action balance equation that accounts for the growth,
propagation, non-linear interaction and dissipation of wind waves in the ocean

by:     

    ∂ N/∂ t + ∇_x · (𝐜_g + 𝐔)N +
    ∂/∂ kk̇N + ∂/∂θθ̇ N = S/σ

where N(k,θ) is the wave action density spectrum, related to the wave
energy density spectrum F(k,θ) where N(k,θ) =
F(k,θ)/σ. Here k is the wavenumber and k̇ its propagation
speed due to depth- or current-induced Doppler shifting, θ is the wave
direction and θ̇ its propagation speed due to depth- or current
refraction, σ is the wave frequency, 𝐜_g is the group
velocity, and 𝐔 the depth-averaged current velocity. On the
right-hand side, S  represents the sum of source terms, including wave growth,
nonlinear interaction and dissipation. WW3 was originally developed on a
regular grid for global operational wave forecasting, with 2-way nesting for
regional applications. More recently, it has been extended to curvilinear grids
for Arctic applications <cit.>, as well as unstructured meshes for
high-resolution coastal application <cit.>. Most
recently, its traditional `card deck' MPI parallel implementation
<cit.> has been supplemented with a more conventional domain-decomposition
 approach using ParaMETIS <cit.> equipped with an optional
 implicit equation solver. Along with these improvements in numeric, source terms suitable for nearshore application have been added, including
depth-induced breaking, reflection, three-wave non-linear interaction, and
wave-ice interaction, amongst others.                    

The WW3 code is written in modular FORTRAN, similar to ADCIRC,
and is broken up into a collection of sub-programs which are run in sequence to
carry out a simulation. The most important of these are ww3_grid, which
compiles the computational mesh, and physics and numerics settings into
mod_def.ww3, a binary resource file, ww3_prep, which preprocesses all
forcing files, ww3_multi, the multi-grid core wave model, and
ww3_ounf and ww3_ounp, which are NetCDF postprocessing routines. To
comply with the ESMF protocol, the core wave model ww3_multi has
been broken up into w3init, w3wave, and w3final to perform the main
steps of model initialization, model advancing, and model finalization,
respectively <cit.>. During the initialization step with
w3init, the configuration of the computational mesh file, including node indices,
geographical location, and the mesh connectivity is read from the binary
resource file mod_def.ww3. Using the domain decomposition from the
PDLIB library, a local and global nodal map is built that references which
nodes reside on which MPI process, and how they are related globally. During the
model advance step w3wave, forcing fields such as water depth, wind
velocity and currents as well as boundary conditions are updated, followed by
the solution of the wave action equation equation by means of fractional
stepping. Output is written to a set of binary output files for later
postprocessing (using ww3_ounf and ww3_ounp). Model finalization is
completed by calling w3final.                       



 §.§ WAVEWATCH III NUOPC cap

The WW3 NUOPC cap carries out the coupling in the three phases
(initialization, advancing and finalizing) described above. A regular grid
version of this cap was developed by <cit.> for global and
regional-scale NUOPC applications. In the present work, this cap was
extended to support unstructured meshes and domain decomposition for
high-resolution coastal modeling. During the initialization of the NUOPC
application, import and export meshes are defined based on the PDLIB
decomposition, a global MPI communicator is created by the ESMF
infrastructure, and a dedicated set of processes are passes to WW3 via a
MPI communicator based on the number of processes requested for WW3 in
the configuration file (Fig. <ref>). During this initialization step, the WW3 cap is
also connected to available import and export field matches accepted by the
communicators.             


After the information exchange among the model components, the
ModelAdvance() subroutine of the WW3 cap calls the w3wave
subroutine to perform the next run interval. Tab. <ref> shows
the list of the exported and imported fields accepted by the WW3 cap for
the current application. It is noted that a larger set of import and export
variables is supported by the WW3 cap in general, including surface
roughness variables, Stokes drift and bed roughness, used in other NUOPC
applications featuring wind waves. For more details, see
<cit.>.



 §.§ Wave-induced stresses

Breaking waves transfer their momentum to ocean currents. Mathematically, this
forcing is expressed in the circulation model as the divergence in the radiation
stresses, as described in some detail below. From spectral wave models such as
WW3, the radiation stress vectors can be evaluated from the computed wave
energy density as follows <cit.>:      

    S_XX   =ρ_w g ∫∫(n-0.5 +ncos^2(θ))F(k,θ)dθ dt 
    
    S_XY   =ρ_w g ∫∫ nsin(θ)cos(θ)F(k,θ)dθ dt 
    
    S_YY   =ρ_w g ∫∫(n-0.5 +nsin^2θ)F(k,θ)dθ dt

where ρ_w is the water density, g is the gravitational acceleration, the
directional wave energy density spectrum F(k,θ) = σ N(k,θ), d
is the water depth, and n is the ratio of the wave group c_g to wave phase
speed c_p for a given depth and frequency, given by:     

n=1/2+kd/sinh(2kd).


In order to account for the impact of the waves on the mean circulation, the
spatial gradient of wave radiation stress per unit area τ_s, waves are
calculated as:

    τ_sx,waves = -( ∂ S_XX/∂ x+∂ S_XY/∂ y)
    τ_sy,waves = -( ∂ S_YY/∂ y+∂ S_XY/∂ x)

and incorporated as additional surface stresses alongside wind stresses and
bottom stresses into ADCIRC's Generalized Wave Continuity Equation and
vertically-integrated momentum equations, following <cit.>.    




§ MODEL SETUP



We utilized the existing Hurricane Surge On-demand Forecast System

unstructured triangular mesh as the base of the
computational domain for the coupled model setup. The HSOFS mesh covers the entire Gulf of Mexico
and extends into the Atlantic Ocean to the approximate longitude of 65∘W, allowing for appropriate
generation of storm surge from atmospheric effects over a large region. The HSOFS 
mesh has 1.8 M nodes and covers the shallow coastal regions up to a 
topographic height of 10 m above local mean sea level with the mesh resolution of approximately 
250 m (See Sec. <ref>.         








    




















 



§ RESULTS



The importance of dynamical coupling of surge and surface waves on the spatial extent
of the inundation and active wave action area were investigated for Hurricane Ike, 2008.

Hurricane Ike was a powerful tropical cyclone that swept through portions of the
Greater Antilles and Northern America in September 2008, wreaking havoc on
infrastructure and agriculture, particularly in Cuba and Texas. The ninth
tropical storm, fifth hurricane, and third major hurricane of the 2008 Atlantic
hurricane season, Ike developed from a tropical wave west of Cape Verde on
September 1 and strengthened to a peak intensity as a Category 4 hurricane over
the open waters of the central Atlantic on September 4 as it tracked westward.
Several fluctuations in strength occurred before Ike made landfall on eastern
Cuba on September 8. The hurricane weakened prior to continuing into the Gulf of
Mexico, but increased its intensity by the time of its final landfall on
Galveston, Texas on September 13.  

The wave-surge coupled application (hearafter Fully coupled) and stand
alone models (hearafter “Stand alone”) were forced with an identical HWRF
meteorological forcing (See Sec. <ref>). 
As a reference, the ADCIRC model results forced with
tidal boundary condition is also presented (heresafter “Only tide”).

The map of the maximum surge level during the whole simulation for the
Fully coupled case is shown in Fig. <ref>a. The maximum surge 
level is calculated by subtracting tidal water level from Fully coupled results.
The results reveal that the most severe inundation during Ike, 2008 with more than 6 m above the
maximum tide level was took place on the east side of the hurricane track in the
region between the Galveston Bay, Tx and the Sabine Lake, Tx. 




The maximum wave contribution to total water level is calculated by subtracting 
Stand alone water level from Fully coupled results (Fig. <ref>b).
It is shown that some of the wave induced inundation (wave setup) occurs at the edge of the
maximum surge where the atmospheric wind setup and negative pressure are at
their maximum strength. This is significant because it shows how wave induced momentum
released from breaking waves increases the total water level which in turn causes the wave active breaking
region to advance further into the land and release wave action on the structures
further landwards. 

Wave height significantly enhanced due to the dynamical coupling of the surge and
wave components. Fig. <ref>a, b represents the maximum wave
height difference during the whole storm between Fully coupled and Stand 
alone (which includes no tidal or surge/inundation effects)
cases. Wave height increases more than 2 m along the track as well as
in the east of the Galveston Bay at the coastal and landfall region. The
comparison at the 6 quick deployed wave gauges shows significant contribution of
the surge on the eastern side of hurricane track in nearshore region (Fig.
<ref>. 

It should be noted that wave setup contribution seen in Fig. 
<ref>b at the Mississippi and Atchafalaya rivers delta
region occurred hours before the actual landfall therefore the maximum surge and 
wave setup in this region did not happen at the same time. However,
it also points to the possibility of experiencing large swells and rip-current
events hours before hurricane makes its actual landfall.

To further analyze this mechanism, we plot changes of the total water level (TWL)
and wave height for the Fully coupled and Stand alone cases at a transect  
shown in Fig. <ref> (The location of the transect is shown in Fig.
<ref>). We also shown the topobathy values along the transect in a
positive upward vertical coordinate system where Z=0 m located at the local mean
sea level (black line in <ref>b). We plotted High Water Mark
observations in the 1 km radius from the transect with red squares.
The TWL and wave height line plots start from the land (kilometer 0) in which ground level is
almost 6 m above the local mean sea level (latitude ∼ 29.80 N) and continues
towards the ocean to 8 m water depth (kilometer 30) in the ocean side (latitude ∼ 29.55 N).
Fully coupled solution for the total water level continuously show enhancement
over the Stand alone results (Fig.<ref>a). Around the shoreline
where the Z∼0 m and farther off-shore, the Stand alone model show greater
total water level in comparison with the Stand alone solution.
This directly relates to wave height and its dissipation shown in Fig.
<ref>b. 

We see that Stand alone model (forced without water surface elevation from surge
component) shows wave height of almost zero right after waves cross the
shoreline (landwards of the Z=0 m; Kilometer∼24). On the other hand, the
wave height from the coupled model show significant wave height of ∼ 2.2 m in the same
region (landwards of the Z=0 m; Kilometer∼24). This pattern is visible
for the wave height evolution even more landwards. This mechanism in
which greaer total water level (∼0.5 m) leads to the potential for more
active local wave generation and propagation and therefore greater
wave dissipation and release wave action which re-ignite enhancement of total
water level in inundated region is presented here.

We also looked at a timeseries of water level at  P shown as a blue triangle in 
the Fig. <ref>. This point is located very close to the maximum 
spatial extent of the inundated region (Fig.
<ref>) which is helpful to examine time variation of the total
water level for Full coupled and Stand alone cases. The ground level is
also plotted by a black dashed-line which shows that P is not wet before the
storm as it is located ∼3 m above the local mean sea level. During the land
fall both Fully coupled and Stand alone cases produced total water level that
inundated the area. Full coupled shows innudation of ∼1 m above the ground
level which is ∼0.5 m above the water level resulted by Stand alone case at
the time of the peak of the storm.
































  
 
 










    


















§ MODEL SYSTEM VALIDATION


We verified the coupled ADCIRC-WW3 application in a step-by-step manner. In the first
step, we verified all the ESMF  intermediate exchange ESMF_state fields before sending
and after receiving by the other model component. Then we performed a basic
verification using a small setup to make sure the coupled ADCIRC-WW3 application runs 
smoothly. Finally, we switched to the full scale inundation test case for
Hurricane Ike, 2008.











 §.§ Laboratory Case (Boers 1996)

To validate the ADCIRC-WW3 coupled system, we first
compared the coupled system results against the laboratory flume
experiment of <cit.>. This experiment was carried out at Delft
University of Technology in Spring 1993 to investigate the interaction 
between wind waves and the mean circulation via wave dissipation
and radiation stress transfer. The geometry of the flume is shown in Fig.
<ref>c. The bathymetry represents an immovable (concrete)
profile of a typical barred beach along the western Dutch coastline. The
<cit.> experiment features three test cases, 1A, 1B and 1C, 
in which 1A and 1B represent violently-breaking, locally-generated wind waves, 
and 1C represents a
mildly-breaking swell. Tab. <ref> shows the wave height and peak
period parameters of the imposed wave conditions at the wave maker, while Fig.
<ref> shows their corresponsing energy spectra.                  





    












Fig. <ref>a,b show the domain decompositions of 
the WW3 and ADCIRC model components, respectively. In both cases the
decomposition was generated using METIS <cit.>, but recall
that in NUOPC/ESMF they do not necessarily need to match - for WW3 the
decomposition was done for 24 cores, while for ADCIRC it was done for 47 cores. The wave component of the coupled
ADCIRC-WW3 was forced with the observed spectra at the upstream wave maker,
and the water level in the flume is initially set to rest. The WW3
component is configured with a time step of 1 s. The relevant wave physics
processes are depth-induced breaking and three-wave nonlinear interaction, using
the source term formulations of, respectively, Battjes and Janssen (1978), with
γ_BJ = 0.80 and α = 1, and <cit.>.
The ADCIRC model component was run at a time step of 1 s, and forced from rest
by only the coupled wave radiation stresses. The NEMS coupling time step
for the two model components is 1 s.                























Fig. <ref> shows that the modeled significant wave heights and
the wave setup produced by coupled ADCIRC-WW3 application are in good
agreement in terms of significant wave heights and the maximum water surface for
all three cases. For the more energetic wave case 1B, wave heights are somewhat
underestimated offshore of the outer bar. This overestimated dissipation appears to result in
the overestimation of the wave setup over this region leading up to the outer
bar. Conversely, in the milder-breaking swell case 1C we observe an
overestimation of wave heights inshore of the outer bar, presumably due to
insufficient depth-induced breaking. This results in a slight underestimation of
nearshore water levels. Nevertheless, these results show good overall skill of
our coupled ADCIRC-WW3 application in capturing the wave height and wave setup
evolution correctly.          







 §.§ Full scale inundation case


Here we present the results of the full scale inundation study for Hurricane
Ike, 2008. To simulate the combined waves and storm surge for Hurricane Ike, 
we utilized the unstructured triangular mesh
used by NOAA's operational Hurricane Surge On-demand Forecast System (HSOFS). We
applied 8 tidal constituents (M2, K1, O1, P1, Q1, N2, S2, K2) at
the model open boundaries (Figs. <ref>).

 




  §.§.§ Atmospheric forcing


The atmospheric forcing for this study is provided by NOAA's Hurricane Weather
Research and Forecasting (HWRF) modeling system, coupled to the MPIPOM ocean model,
and empowered by a movable multi-level nesting technology
<cit.>. The model grid is triple-nested using telescopic,
two-way interactive horizontal grid resolutions from synoptic with 
0.18^∘ resolution as the outer domain (spanning about 75^∘
× 75^∘),  to moving storm nest with 0.06^∘ resolution
(10^∘ × 10^∘) and core of about 6^∘ ×
6^∘ with 0.02^∘ resolution. These nests follow the hurricane
best track, ensuring the highest resolution around the eye of a hurricane. In
this study, we have interpolated the hourly HWRF model outputs from
multiple cycles initiated with analysis data and 9 forecast time steps. Every 6
hours, reanalysis data from the next cycle are smoothly ramped into the wind and
pressure fields. The atmospheric forcing has been validated against National
Data Buoy Center (NDBC) and satellite altimeter data. We extracted wind velocity
at 10 m height and surface pressure from the original GRIB2 output files
and saved them in NetCDF format. The ATMesh NUOPC data cap reads
the meteorological forcing from NetCDF file and provides it to ADCIRC
and WW3 caps at every coupling time step. The HWRF model was forced
with initial and boundary conditions provided by NOAA's Global Forecast System
(GFS) with 0.5^∘ spatial grid resolution.                      

Fig. <ref> gives an impression of the quality of the HWRF
atmospheric forcing fields used to force the coupled ADCIRC-WW3 model, by
comparing 10 m wind speed and direction against NDBC buoy observations in the
Gulf of Mexico (see Fig. <ref> for locations). We can see that the agreement is generally good, in particular
for wind directions. However, a tendency to overestimate the wind speeds at the
storm peak is found at the mid-Gulf NDBC buoys 42001 and 42002, and the shelf
buoy 42019. The high bias in the wind speed is particularly evident at landfall,
as seen at NDBC 42035, located on the shelf just offshore of Galveston. It would
be expected that these overestimated winds would lead to a degree of
overestimation of the locally-generated significant wave heights. It
is also worth mentioning that the NOAA buoy 42035 broke free during the storm
(See Blockedhttps://www.ndbc.noaa.gov/hurricanes/2008/ike/). Therefore, the
observation time series might not be at the same location of extraction 
coordinate from the models. 
















































 












  §.§.§ ADCIRC model component validation

We validated the total water level results against two sets of water level
observations. The tidal gauge time-series were measured by the NOAA's Center for
Operational Oceanographic Products and Services (CO-OPS) and High Water Marks
(HWM) were measured and provided by United States Geographical Survey (USGS).

Comparison between the total water surface elevation from the coupled ADCIRC-WW3 application and tidal
gauges time series for four stations close to the actual hurricane track (location 
of the tidal gauges in Fig. <ref> are presented in
Fig. <ref>a-d).
The modeled water surface elevation agrees well with the observation in terms of the 
level and timing of the peak of the storm inundation. 





High Water Marks are an important source of observations for validation and
enhancement of the storm surge and flood inundation studies. After significant
flooding due to a land-falling hurricane, a rapid high water mark (HWM) data
collection by USGS takes place to document the event and to help improving
future disaster preparedness activities. Comparison of the total water level from the 
coupled ADCIRC-WW3 application and HWMs observation reveal that both are 
in general agreement in particular around the hurricane track in the land-fall 
region (Fig. <ref>). 





We also compared scatter plots and a number of statistical metrics of the total
water levels and HWM data (Fig. <ref>). Both Fully coupled and 
Stand alone results show
underestimation of the model in comparison to observations. However, upon 
inclusion of the wave forcing, an improvement in the error statistics is found, such 
that a reduction in the relative bias (RB) from -0.594 m to -0.392 m, and a reduction 
in root mean square error (RMSE) from 0.899 m to 0.832 m were presented.




 


 

 
 




  §.§.§ WW3 component verification

The coupled ADCIRC-WW3 coupled application performance with respect to the
wave modeling aspects was firstly assessed in the offshore and over the shelf of the Gulf of Mexico using 
NDBC monitoring buoys, and secondly at the location of landfall by 
using rapid-deployment pressure sensors, placed by <cit.>.

Considering the NDBC monitoring buoys, four groups can be distinguished. At mid-Gulf 
(NDBC 42001 and 42002), we find some overestimation of significant wave heights 
during the passing of the hurricane in response to the somewhat overestimated 
wind speeds seen above (Fig. <ref>a). Similarly, the peak wave periods 
are overestimated, reflecting excessive wind-sea growth (Fig. <ref>b). 
However, wave direction is reproduced well (Fig. <ref>c). Interestingly, 
in the far-field, towards the eastern half of the Gulf (NDBC 42036, 42039 42040), 
significant wave heights tended to be underestimated during the passing of the 
storm, while the peak period and direction were captured well. Moving onto the 
shelf in the region of landfall (NDBC 42019 and 42020) the agreement of the model 
results with the observations improved in general, although there is still 
overestimation in the significant wave height and peak period. Finally, at 
NDBC 42035, just offshore of the landfall location at Galveston, all wave model 
parameters are in good agreement with the observations at the storm peak. However, 
the significant wave height is underestimated at this shallow water location 
just ahead of the storm peak which could be related to the omission of the
forerunner effect in the coupled surge model component. Regarding the differences between the 
coupled and stand-alone versions of the wave model, these can only be seen, 
as expected, at the two shallower stations NDBC 42020 and 42035. As expected, 
the increased water levels seen above increase the modeled 
significant wave height at these stations.

It should be noted that our concentration here is to ascertain general
performance of the coupled model. In terms of forerunner, we reference our reader to
<cit.>. We diagnose the forerunner surge as being generated 
by Ekman setup on the wide and shallow shelf. The longer forerunner time 
scale additionally served to increase water levels significantly in 
narrow-entranced coastal bays. 






  


Moving more towards shoreline, model results at the pressure sensors deployed 
in the surf zone by <cit.> show a significant sensitivity
 to the inclusion of variable water levels from the coupled ADCIRC model (Fig.
 <ref>). Referring to Fig. <ref>b for the station
 locations, the most significant of these are the stations ANDKNDY-X, Y and Z
 located under the eastern half of the landfalling hurricane. For these three
 stations, we can see a large influence of the added surge level on the
 significant wave height at the storm peak, in all cases improving the agreement
 with observations. 
 We note that, as discussed above, the stations U to Z also might 
 relate to underestimation of the water surface elevation during forerunner.
 Stations S and R, located to the
 south of the landfall location received mostly offshore winds and a water level
 set-down. The significant wave heights at these stations are reproduced well.          

  














Fig. <ref> provide a more in-depth view of the results at 
the <cit.> nearshore stations by considering their spectrograms. 
This figure shows a number of important features of the nearshore wave field 
transformation under both the coupled and stand-alone models. Up to September 
13, the simulated variance density at the northerly stations ANDKNDY-X, Y and Z 
is underestimated by both the stand-alone (Fig. <ref>b) and coupled 
(Fig. <ref>c) models in connection with non-resolved forerunner
effects.




After September 13, the variance density abruptly increases, as the main storm peak 
arrives with landfall. At this time, the coupled model (panel c) shows greater 
levels of variance density than the stand-alone model (panel b), which agrees 
better with the observations. A final important feature of the nearshore spectra 
is the frequency upshift of the spectral peak, due to the nonlinear 
three-wave interactions. This strong upshift at stations V, W, X, Y, Z is not 
seen as intensely in the model results, indicating an underestimation 
of the magnitude of this process by the LTA nonlinear interaction source term. 
Interestingly, the stand-alone model captured this upshift process somewhat 
better than the coupled model, since the depth underestimation in the former 
fortuitously enhances the computed nonlinear interaction (see stations X, Y, Z).

 
  
 


















§ SUMMARY AND CONCLUSIONS



We developed a flexible coupling application for coastal inundation studies in the framework of
NEMS using NUOPC/ESMF infrastructure. This application includes
NUOPC model interfaces (or caps) for ADCIRC and WW3, and also a
data cap to read and provide atmospheric forcing fields from HWRF model
results. The application is examined using a standard laboratory flume case for 
set-up and wave dissipation, and validated for Hurricane Ike, 2008 on NOAA's
HSOFS mesh, a 1.8 M node triangular mesh with a nominal resolution of ∼500
m. The model skills and improvement due to wave effects on the final inundation were examined and discussed using
time series from tide gauges and high water marks observations.

In general, it can be concluded that the coupling improves the performance of
the ADCIRC and WW3 model components. The <cit.> case
shows that the coupling between these two models behaves correctly under laboratory conditions. In the Ike field case, the simulated water
levels of ADCIRC showed generally better agreement with observations upon
inclusion of the coupled wave effects from WW3. Accounting for the
forerunner effect is expected to improve the results further. Conversely, the 
coastal wave field simulated by WW3 improves significantly
upon coupling with ADCIRC, in particular in the surf zone and in inundated
regions.

As outlined in the Introduction, the first application of the NUOPC model
combination (or App) described here is the Named Storm Event Model (NSEM), a high-resolution hindcast model
which is being developed to meet the requirements of the COASTAL Act (2012). However, we
anticipate that these flexible and generic NUOPC coupling caps for ADCIRC and WW3 will
enable future development and seamless inclusion of various additional model components for
forecasting applications. For example, the current cap development provides the possibility to
perform three-way dynamical coupling of ADCIRC, WW3 and atmospheric
prediction models.
In the hindcast application described here, we limited the atmospheric model, run offline, to a
one-way forcing to ADCIRC-WW3 via a data cap. However, the ADCIRC and
WW3 caps are capable of dynamically importing atmospheric model forcing,
and exporting current velocity, water surface elevation (accounting for inundated zones which alter the surface roughness), and
enhanced sea surface roughness due to the presence of waves (Charnock parameter) back to the
connected atmospheric model component.

Further future coastal applications of the presented flexible NUOPC
modeling framework are incorporating processes such as river and inland flooding and sea ice. As part of the Named
Storm Event Model, work is currently under way to incorporate river flow into this coupled
system by including a NUOPC cap for WRFHydro. This
extension will allow the modeling of complex coastal flooding events such as Hurricane Harvey (2017), which featured
extreme rainfall and inland flooding alongside the surge and wave forcing from the ocean. Sea
ice has significant and complex impacts on coastal surge
<cit.> as well as dissipative and dispersive effects on the coastal 
wave field <cit.>. The
flexible coastal modeling system presented here can be extended to include these coupled processes via
an existing NUOPC cap for the sea ice model CICE
<cit.>.




§ ACKNOWLEDGEMENTS

This work has been supported by The Consumer Option for an Alternative System to Allocate Losses
(COASTAL) Act Program and Water Initiative project within the National
Oceanic and Atmospheric Administration (NOAA). The WaveWatchIII development has been supported by the US. 
Army Corps of Engineers (USACE). The authors acknowledge Dr. Andrew Kennedy for providing the
in-situ data. The results presented here are preliminary and for presenting the
incremental developments of coupled application which is under development to
serve above mentioned projects. The authors also would like to acknowledge Dr.
Chris Massey for his constructive communications.



spbasic



25
natexlab#1#1









[#1],#1
[Ardhuin and Roland(2013)]ardhuin2013development
authorArdhuin, F., authorRoland, A.,
  year2013.
titleThe development of spectral wave models: coastal and
  coupled aspects, in: booktitleProceedings of Coastal Dynamics,
  pp. pages25–38.

[Biswas K et al.(2018)Biswas K, Bernardet, Abarca, Ginis, Grell,
  Kalina, Kwon, Liu, Liu, Marchok, Mehra, Newman, Sheinin, Sippel, Subramanian,
  Tallapragada, Thomas, Tong, Trahan, Wang, Yablonsky, Zhang and
  Zhang]hwrf2017
authorBiswas K, M., authorBernardet, L.,
  authorAbarca, S., authorGinis, I.,
  authorGrell, E., authorKalina, E.,
  authorKwon, Y., authorLiu, B., authorLiu,
  Q., authorMarchok, T., authorMehra, A.,
  authorNewman, K., authorSheinin, D.,
  authorSippel, J., authorSubramanian, S.,
  authorTallapragada, V., authorThomas, B.,
  authorTong, M., authorTrahan, S.,
  authorWang, W., authorYablonsky, R.,
  authorZhang, X., authorZhang, Z.,
  year2018.
titleHurricane Weather and Research and Forecasting (HWRF)
  Model: 2017 scientific documentation.
typeTechnical Report. Development Testbed Center, UCAR.

[Boer(1996)]boer1996Surfzone
authorBoer, M., year1996.
titleSimulation of a surf zone with a barred beach. Report
  1. Wave heights and wave breaking.
typeTechnical Report. Dept. of Civil Engineering, Delft
  University of Technology.

[Campbell and Whitcomb(2013)]campbell2013coupling
authorCampbell, T.J., authorWhitcomb, T.R.,
  year2013.
titleCoupling Infrastructure and Interoperability Layer
  Extension Across All Earth System Components.
typeTechnical Report. NAVAL RESEARCH LAB STENNIS
  DETACHMENT STENNIS SPACE CENTER MS.

[Dietrich et al.(2011)Dietrich, Zijlema, Westerink, Holthuijsen,
  Dawson, Luettich, Jensen, Smith, Stelling and Stone]dietrich2011modeling
authorDietrich, J., authorZijlema, M.,
  authorWesterink, J., authorHolthuijsen, L.,
  authorDawson, C., authorLuettich, R.A.,
  authorJensen, R., authorSmith, J.,
  authorStelling, G., authorStone, G.,
  year2011.
titleModeling hurricane waves and storm surge using
  integrally-coupled, scalable computations.
journalCoastal Engineering volume58,
  pages45–65.

[Eldeberky and Battjes(1996)]eldeberky1996spectral
authorEldeberky, Y., authorBattjes, J.A.,
  year1996.
titleSpectral modeling of wave breaking: Application to
  boussinesq equations.
journalJournal of Geophysical Research: Oceans
  volume101, pages1253–1264.

[Gochis et al.(2013)Gochis, Yu and Yates]gochis2013wrf
authorGochis, D., authorYu, W., authorYates,
  D., year2013.
titleThe wrf-hydro model technical description and
  user’s guide, version 1.0.
journalNCAR technical document volume120.

[Hill et al.(2004)Hill, DeLuca, Balaji, Suarez and
  Silva]hill2004architecture
authorHill, C., authorDeLuca, C.,
  authorBalaji, authorSuarez, M.,
  authorSilva, A.d., year2004.
titleThe architecture of the earth system modeling
  framework.
journalComputing in Science & Engineering
  volume6, pages18–28.

[Jacob et al.(2005)Jacob, Larson and Ong]jacob2005m
authorJacob, R., authorLarson, J.,
  authorOng, E., year2005.
titleM× N communication and parallel interpolation
  in Community Climate System Model version 3 using the Model Coupling
  Toolkit.
journalInternational Journal of High Performance Computing
  Applications volume19, pages293–307.

[Karypis(2011)]karypis2011metis
authorKarypis, G., year2011.
titleMETIS and ParMETIS, in:
  booktitleEncyclopedia of parallel computing.
  publisherSpringer, pp. pages1117–1124.


[Kennedy et al.(2011)Kennedy, Gravois, Zachry, Westerink, Hope,
  Dietrich, Powell, Cox, Luettich and Dean]kennedy2011origin
authorKennedy, A.B., authorGravois, U.,
  authorZachry, B.C., authorWesterink, J.J.,
  authorHope, M.E., authorDietrich, J.C.,
  authorPowell, M.D., authorCox, A.T.,
  authorLuettich, R.A., authorDean, R.G.,
  year2011.
titleOrigin of the hurricane ike forerunner surge.
journalGeophysical Research Letters volume38.

[Lemmen et al.(2017)Lemmen, Hofmeister, Klingbeil, Nasermoaddeli,
  Kerimoglu, Burchard, Kösters and Wirtz]lemmen2017modular
authorLemmen, C., authorHofmeister, R.,
  authorKlingbeil, K., authorNasermoaddeli, M.H.,
  authorKerimoglu, O., authorBurchard, H.,
  authorKösters, F., authorWirtz, K.W.,
  year2017.
titleModular system for shelves and coasts (mossco v1.
  0)-a flexible and multi-component framework for coupled coastal ocean
  ecosystem modelling.
journalarXiv preprint arXiv:1706.04224 .

[Luettich Jr et al.(1992)Luettich Jr, Westerink and
  Scheffner]luettich1992adcirc
authorLuettich Jr, R.A., authorWesterink, J.J.,
  authorScheffner, N.W., year1992.
titleADCIRC: An Advanced Three-Dimensional Circulation
  Model for Shelves, Coasts, and Estuaries. Report 1. Theory and Methodology of
  ADCIRC-2DDI and ADCIRC-3DL.
typeTechnical Report. COASTAL ENGINEERING RESEARCH CENTER
  VICKSBURG MS.

[Moghimi et al.(2012)Moghimi, klingbeil, Graewe and
  Burchard]moghimi2012three
authorMoghimi, S., authorklingbeil, K.,
  authorGraewe, U., authorBurchard, H.,
  year2012.
titleThree dimensional wave current interaction in
  nearshore regions, in: booktitleJONSMOD2012 , Brest, France.

[Moghimi et al.(2019)Moghimi, Vinogradov, Myers, Funakoshi, Van der
  Westhuysen, Abdolali, Ma, Liu et al.]moghimi2019development
authorMoghimi, S., authorVinogradov, S.,
  authorMyers, E.P., authorFunakoshi, Y.,
  authorVan der Westhuysen, A.J., authorAbdolali, A.,
  authorMa, Z., authorLiu, F., et al.,
  year2019.
titleDevelopment of a Flexible Coupling Interface for
  ADCIRC Model for Coastal Inundation Studies.
typeTechnical Report. NOAA National Ocean Service.

[Rogers et al.(2018)Rogers, Posey, Li and
  Allard]rogers2018forecasting
authorRogers, W., authorPosey, P., authorLi,
  L., authorAllard, R., year2018.
titleForecasting and Hindcasting Waves In and Near the
  Marginal Ice Zone: Wave Modeling and the ONR Sea State Field Experiment.
typeTechnical Report. NAVAL RESEARCH LAB STENNIS
  DETACHMENT STENNIS SPACE CENTER MS STENNIS SPACE ….

[Roland et al.(2012)Roland, Zhang, Wang, Meng, Teng, Maderich,
  Brovchenko, Dutour-Sikiric and Zanke]roland2012fully
authorRoland, A., authorZhang, Y.J.,
  authorWang, H.V., authorMeng, Y.,
  authorTeng, Y.C., authorMaderich, V.,
  authorBrovchenko, I., authorDutour-Sikiric, M.,
  authorZanke, U., year2012.
titleA fully coupled 3d wave-current interaction model on
  unstructured grids.
journalJournal of Geophysical Research: Oceans
  volume117.

[Tallapragada et al.(2014)Tallapragada, Kieu, Kwon, Trahan, Liu, Zhang
  and Kwon]tallapragada2014evaluation
authorTallapragada, V., authorKieu, C.,
  authorKwon, Y., authorTrahan, S.,
  authorLiu, Q., authorZhang, Z.,
  authorKwon, I.H., year2014.
titleEvaluation of storm structure from the operational
  hwrf during 2012 implementation.
journalMonthly Weather Review volume142,
  pages4308–4325.

[Theurich et al.(2016)Theurich, DeLuca, Campbell, Liu, Saint,
  Vertenstein, Chen, Oehmke, Doyle, Whitcomb et al.]theurich2016earth
authorTheurich, G., authorDeLuca, C.,
  authorCampbell, T., authorLiu, F.,
  authorSaint, K., authorVertenstein, M.,
  authorChen, J., authorOehmke, R.,
  authorDoyle, J., authorWhitcomb, T., et al.,
  year2016.
titleThe earth system prediction suite: toward a
  coordinated us modeling capability.
journalBulletin of the American Meteorological Society
  volume97, pages1229–1247.

[Tolman(2002)]TOLMAN200235
authorTolman, H.L., year2002.
titleDistributed-memory concepts in the wave model
  WAVEWATCH III.
journalParallel Computing volume28,
  pages35 – 52.




[Tolman et al.(2009)]tolman2009user
authorTolman, H.L., et al., year2009.
titleUser manual and system documentation of wavewatch iii
  tm version 3.14.
journalTechnical note, MMAB Contribution
  volume276, pages220.

[Valcke et al.(2012)Valcke, Balaji, Craig, DeLuca, Dunlap, Ford,
  Jacob, Larson, O'Kuinghttons, Riley et al.]valcke2012coupling
authorValcke, S., authorBalaji, V.,
  authorCraig, A., authorDeLuca, C.,
  authorDunlap, R., authorFord, R.,
  authorJacob, R., authorLarson, J.,
  authorO'Kuinghttons, R., authorRiley, G., et al.,
  year2012.
titleCoupling technologies for earth system modelling.
journalGeoscientific Model Development
  volume5, pages1589.

[Warner et al.(2008)Warner, Sherwood, Signell, Harris and
  Arango]warner2008development
authorWarner, J., authorSherwood, C.,
  authorSignell, R., authorHarris, C.,
  authorArango, H., year2008.
titleDevelopment of a three-dimensional, regional,
  coupled wave, current, and sediment-transport model.
journalComputers & Geosciences volume34,
  pages1284–1306.

[Westerink et al.(2018)Westerink, Joyce, Grumbine, van der Westhuysen,
  Feyen, Pringle and Wirasaet]westerink2018coupled
authorWesterink, J.J., authorJoyce, B.,
  authorGrumbine, R., authorvan der Westhuysen, A.,
  authorFeyen, J., authorPringle, W.J.,
  authorWirasaet, D., year2018.
titleCoupled tides, storm surge and waves under varying
  ice coverages along alaska’s bering and chukchi coasts.
journalCoastal Engineering Proceedings
  volume1, pages69.

[Zhang et al.(2016)Zhang, Gopalakrishnan, Trahan, Quirino, Liu, Zhang,
  Alaka and Tallapragada]zhang2016representing
authorZhang, X., authorGopalakrishnan, S.G.,
  authorTrahan, S., authorQuirino, T.S.,
  authorLiu, Q., authorZhang, Z.,
  authorAlaka, G., authorTallapragada, V.,
  year2016.
titleRepresenting multiple scales in the hurricane weather
  research and forecasting modeling system: Design of multiple sets of movable
  multilevel nesting and the basin-scale hwrf forecast application.
journalWeather and Forecasting volume31,
  pages2019–2034.





§ APPENDIX A: METRICS FOR THE EVALUATION OF DATA-MODEL AGREEMENT



In order to assess model performance for water levels, root mean square error
(RMSE), BIAS, relative BIAS (RB), Correlation (Cor), Index of Agreement (IA) and
peak error (Peak) were used.

The RMSE is given by

    RMSE = √(1/N∑_i=1^N (M_i - O_i)^2)

where M_i is the modeled data, O_i is the measured data and N is the total number of data.

BIAS shows the systematic deviation from the observations and is given by

    BIAS = ∑_i=1^N (M_i - O_i)


Relative BIAS (RB) shows relative systematic deviation from the observations
and is given by

    RB = ∑_i=1^N (M_i - O_i)/N ⟨ O_i ⟩



Peak error is calculated by

    PEAK = maxO - maxM



The Index of Agreement (IA) is formulated as

    IA = 1 - ∑_i=1^N (M_i - O_i)^2/∑_i=1^N ((| M_i - ⟨ O ⟩| + | O_i - ⟨ O ⟩|) ^2 )

where brackets, ⟨·⟩, denote time averaging. IA=1 shows perfect agreement and IA=0 means complete disagreement.

The Pearson correlation (Cor) coefficient is calculated by

    Cor =∑_i=1^N (O_i -  ⟨ O ⟩) (M_i -  ⟨ M ⟩) /√(∑_i=1^N (O_i -  ⟨ O ⟩)^2  ∑_i=1^N
     (M_i -  ⟨ M ⟩)^2 )

It has a value between +1 and −1, where 1 is total positive linear
correlation, 0 is no linear correlation, and −1 is total negative linear correlation 

