Development and validation of a coastal ocean forecasting system for Puerto Rico and the U.S. virgin islands

M. Solno , M. Cnls , S. Leonrdi , ?

a Naval Research Laboratory Stennis Space Center, MS 39529-5004

b Center for Ocean Science and Engineering, Department of engineering Sciences and Materials,Caribbean Coastal Ocean Observing System, University of Puerto Rico at Mayaguez

c Dept. Mechanical Engineering, The University of Texas at Dallas, TX USA

Abstract An ocean forecasting system has been developed for the coastal area of Puerto Rico and the U.S. Virgin Islands. This paper presents the development and validation of the Puerto Rico Operational Regional Ocean Modeling System (PRO-ROMS) developed by the Caribbean Coastal Ocean Observing System (CARICOOS), which aims at improving currently available short-term forecasts of ocean currents for this region, and to resolve sub-mesoscale variability absent from the currently operational regional systems. The coastal ocean forecasting system is based on the Regional Ocean Modeling System (ROMS), and provides a 3-day forecast of ocean currents, water levels, temperature and salinity. Initial and lateral boundary conditions are derived from the U.S. Naval Oceanographic Office (NAVOCEANO) operational AmSeas forecasting system, a 3.6-km resolution application of the regional Navy Coastal Ocean Model (NCOM) that encompasses the Gulf of Mexico and Caribbean Sea. Meteorological conditions are interpolated from the Navy’s Coupled Ocean-Atmosphere Mesoscale Prediction System(COAMPS) with the exception of surface stresses, which are computed from a 2-km application of the Weather Research and Forecasting(WRF) model used by National Centers for Environmental Protection’s (NCEP) National Digital Forecast Database (NDFD). Tidal variability is imposed using ROMS spectral forcing by specifying the harmonic phases and amplitudes from the TPXO global inverse tide solutions at the boundaries and sub-tidal conditions are imposed by filtering out high frequencies of barotropic variables from the lateral boundary conditions interpolated from AmSeas. Modeled results of sea surface height are validated with NOAA tide gauges, ocean currents are validated with Acoustic Doppler Current Profilers (ADCP) and High Frequency Radars (HFR). Computed statistics show improvement of the ocean current forecast in PRO-ROMS compared to AmSeas, especially at higher frequencies dominated by the tides and in small channels where the bathymetry is better represented by the higher resolution coastal ocean forecasting system.

Keywords: Ocean modeling; Operational oceanography; Forecasting.

1.Introduction

Ocean modeling and forecasting has become an essential tool used by regional and international ocean observing systems in an attempt to provide an accurate estimate of the ocean state. The continuous decrease in the cost of computational resources has allowed regional associations, dedicated to ocean observing, to develop high-resolution coastal ocean models which aim at improving forecasts of existing operational systems for the domain of interest. This paper presents the development and validation of the Puerto Rico Operational Regional Ocean Modeling System (PRO-ROMS), an ocean forecasting system for the coastal region of Puerto Rico and the U.S. Virgin Islands. This system has been developed for the Caribbean Coastal Ocean Observing System (CARICOOS), one of eleven coastal observing systems and regional associations which along with federal agencies constitute the national coastal component of the U.S. Integrated Ocean Observing System (IOOS).

At present, ocean circulation models available for the coastal areas of Puerto Rico provide forecasts with horizontal resolutions up to 1/30 degrees (3.6km) and 3-hour output frequency. Most of the previous ocean circulation studies in Puerto Rico have focused on barotropic tides (e.g., Kjervfe[1] , Torres and Tsimplis [2] ) or water transport through island passages(e.g., Wilson [3] , Johns [4] ), but provide little information on the three-dimensional structure of the ocean currents. The generation of high-amplitude coastal seiches in Puerto Rico was previously investigated by Giese and Chapman [5] , in which they confirm tidal forcing as a necessary requisite for their formation and further argue that the seiches at this site are excited by deep-sea tide-generated internal waves travelling from the Aves ridge, and later confirmed by MODIS Terra/Aqua satellite observations (Giese 1989). In a recent study performed by Gonzalez (2015) [6] , a two-dimensional high-resolution (O(100 m)) ocean circulation model was developed and validated with good agreement in sea surface height (SSH) when compared to observations at NOAA tide gauges in the near-shore coastal shelf region, and suggests that local wind forcing may play a role on the formation of seiches. A study of the Caribbean Current variability by Richardson [7] using satellite-tracked surface drifters identify the presence of two bands of jets in the eastern Caribbean travelling westward. The formation of cyclonic eddies near passages through the eastern islands of the Antilles and others emerging inside the islands suggests that these cyclones are often formed by flows through passages and may be intensified by local wind circulation.

The aim of the coastal model presented here is two-fold: to simulate higher frequency phenomena present in coastal shelf regions, providing a database of high-resolution 3D ocean currents, temperature and salinity for use by the scientific community and to improve the short term forecast of ocean currents in shallow shelf regions, enhancing observing and forecasting capabilities for CARICOOS and the community of Puerto Rico and the U.S. Virgin islands. To achieve this,as normally done in coastal models, we use a downscaling methodology to propagate the large scale dynamics obtained from a coarser, large scale model through the boundaries in combination with a high resolution grid to simulate the higher frequency dynamics present in coastal areas. The use of finer grids implies that smaller space and time scales are resolved,but because the open boundary condition problem is ill-posed,errors at the boundary may propagate inside yielding no improvement over the parent model.

The performance of ocean models solving the primitive equations in a limited domain is significantly constrained by the specification of the open boundaries and model forcing. In configuring open boundary conditions for realistic regional applications using split-explicit ocean models such as ROMS, some considerations should be given to dynamical influences which change considerably depending on the underlying physical processes characteristic to the region under study, most importantly those occurring at the open boundaries. One of these dynamical influences is the dominant scale of ocean variability, namely meso and submesoscales. Submesocale currents (SMC) have an approximate scale range of 0.1–10 km in the horizontal direction, 0.01–1 km in the vertical direction and time scales of hours to several days, which are the typical space-time scales in coastal ocean forecasting systems. SMC emerge spontaneously in coastal ocean models when the horizontal grid resolution is ofO(1 km )(McWilliams, 2017 [8] ). From a physical standpoint, the transition from mesoscale to submesoscale typically occurs atRo＞ 1, when geostrophic and hydrostatic-momentum balance start to break down. Submesoscale currents are generated through instabilities, especially in the vicinity of surface and bottom boundary layers (McWilliams 2017 [8] ), ubiquitous in shelf seas.

Although there seems to be no clear transition from mesoscale to submesoscale (Capet et al. 2008 a,b,c [9–11] )and coastal phenomena, some important distinctions can be made between coastal and regional ocean models. Firstly,coastal ocean models deal with the transition from the deep ocean (depths larger than 1000 m) to shelf seas (depths between 100–10 m), leading to significant energy dissipation via bottom friction and internal tide generation and breaking. The interaction of ocean currents and topography affects the behavior of surface elevation, which may propagate in the form of trapped waves or other higher frequency phenomena difficult to elucidate in coarser regional models. Secondly, tideinduced currents are more energetic in shallow waters than in the deep ocean, which at similar depths result in higher frequency variability when compared to the open ocean (Haidvogel, 2000 [12] ).

The most influential factor in short term forecasts (3-day) is the model’s initialization. Global and regional ocean forecasting systems usually employ some data assimilation method in which several data streams such as satellite andinsitusea surface temperature, satellite altimeter-derived sea level anomaly, temperature and salinity profiles among others are used to iteratively improve short-range forecasts. However, data assimilation algorithms rely on the assumption that the background field is already a good approximation,otherwise the system will perform large corrections to the model which may yield detrimental results. In theory, higherresolution models provide improved dynamical solutions as sub-grid scales decrease, resolving smaller eddies and less parameterization of turbulent processes are needed. In practice, higher-resolution systems will provide analyses with enhanced spatial detail but will be less skillful at predicting the evolution of mesoscale eddies (Sandery et al., [13] ).

Implementation of data assimilation methods in coastal zones becomes more complicated, as observations are usually more sparse due to the much smaller space and time scales found in coastal areas, which is problematic for the computation of the error covariance matrix used for the interpolation and smoothing of sparse observations (Kourafalou et al.,2015 [14] ). Consequently, implementation of data assimilation methods for high resolution coastal applications requires the availability of observations sampled at similar time and space scales to that of the model, as well as previous implementation and validation of the coastal model, to ensure that model dynamics, initial and boundary conditions are configured in such a way that the background field is the best possible estimate of the ocean state. The present implementation of the operational model does not assimilate measurements, but lays the foundation for such an implementation.

The following section provides a brief description of the ROMS model with a detailed explanation of the downscaling methodology, including boundary conditions, model resolution, bathymetry smoothing, meteorological and tidal forcing.Then the validation of the modeled results with observations is presented in Section 3 , where simulated water levels and ocean currents are compared with hydrographic measurements at the shelf and shelf edge. Model performance is summarized using the Taylor diagram [15] , which condenses centered root-mean square error, correlation coefficient, forecast and observed variances into one diagram for a quantitative comparison. A qualitative assessment of the spatial distribution of ocean currents is performed by comparing the modeled surface currents with High-Frequency-Radar (HFR) observations. A grid sensitivity analysis is performed to assess the forecast skill of the model by increasing the uniform horizontal resolution of the operational ROMS grid is changed using a refinement coefficient over the parent NCOM model of 1, 3 and 5, corresponding to 3.6 km, 1.2 km and 0.72 km grid resolutions respectively, and results are briefly reported in section 3.5.. The last subsection of the results presents a particle tracking experiment to provide further validation of PRO-ROMS in simulating and predicting the trajectories of passive surface tracers. Finally, significant results are summarized in the last section, emphasizing overall model performance, improvements over the AmSeas regional model as well as constraints and limitations of the coastal model.

2.Model configuration

2.1. Governing equations

The operational system is based on the Regional Ocean Modelling System (ROMS), a primitive equation ocean model widely used by the oceanographic community to simulate and predict ocean currents. ROMS solves the Reynolds Averaged Navier-Stokes equations, and uses a split-explicit time step to advance the depth-averaged shallow water equations with a short time step, and a longer time step is used for the baroclinic mode (3D) for computational efficiency. The Boussinesq approximation is assumed such that density variations in the horizontal direction are neglected and under the hydrostatic approximation the buoyancy force balances the vertical pressure gradient. For a detailed explanation of the numerical code the reader is referred to Shchepetkin and McWilliams(2005) [16] .

The advective terms in the momentum and passive tracer equations are solved using a 3rd order upstream predictorcorrector scheme for the horizontal direction (Shchepetkin and McWilliams, 1998 [17] ), which allows the generation of subscale geostrophic turbulence in high-resolution coastal models. For the vertical direction a 4th order centered scheme is used with vertical interpolation by conservative splines. The vertical turbulence uses the General Length Scale framework,with ak?∈scheme to solve the closure problem. Under the GLS framework, a smoothing function of buoyancy/shear is applied in the horizontal equation and a polynomial smoothing is applied to the computation of the bulk Richardson number.

2.2. Domain and bathymerty

The computational domain is approximately 400 km ×200 km and encompasses the shelf region of the island of Puerto Rico and the U.S. Virgin Islands, located in the middle of the Antilles Archipelago. The model bathymetry, shown in Fig. 1 , is generated from the SRTM30 PLUS database (Becker et al., 2009 [18] ), which combines a wide variety of data sources to produce an accurate representation of the ocean floor at 30 second-arc resolution, which is roughly the same as the model’s horizontal resolution. The northern boundary intersects the Puerto Rico trench which extends for 800 km longitudinally and surpasses depths of 8000 m; however the model bathymetry is truncated at 5000 m to avoid extrapolations when interpolating lateral boundary conditions from the parent model.

2.3. Grid

A rectangular orthogonal grid with uniform 1/90 degree(1.2 km) resolution is used for the horizontal direction, and 32 non-uniform levels are used in the vertical direction. The vertical grid is stretched towards the surface and slightly towards the bottom to better capture the top and bottom boundary layers. The critical depth ( hc) is set to 100 m, which is approximately the depth of the shelf and the pycnocline in the region (Hernandez-Guerra and Joyce, 2000 [19] ; Hu et al.,2004 [20] ), and determines the depth at which the vertical grid becomes uniform in order to solve the top and bottom boundary layer interaction.

2.4. Downscaling and initialization

Currently, the system is initialized daily from the U.S.Naval Oceanographic Office (NAVOCEANO) operational AmSeas model forecast, a 1/30 degree resolution application of the regional Navy Coastal Ocean Model (NCOM) that encompasses the Gulf of Mexico and Caribbean Sea. The regional NCOM ocean prediction system assimilates all qualitycontrolled observations in the region including satellite sea surface temperature and altimetry, as well as surface and profile temperature and salinity data using the Navy Coupled Ocean Data Assimilation (NCODA) system. Boundary conditions are applied to the AmSeas model from the NAVOCEANO operational 1/12 degree Global HYCOM. The Am-Seas system provides a 4-day forecast of water levels, ocean currents, temperature and salinity at 3-hour intervals.

Fig. 1. SRTM30 PLUS model bathymetry with 30 arc-second horizontal resolution (approx. 1 km). Isobaths are specified for ?10 m, ?100 m, ?500 m,?1000 m, ?2000 m, ?3000 m, ?4000 m and ?5000 m to show the steep slopes at the shelf edges. Land areas are masked in white.

Downscaling from the regional NCOM model is performed using a single one-way nesting approach designed to minimize the child-parent error at the boundaries and optimized to capture high-frequency signals important in coastal circulation. Initial and lateral boundary conditions for the baroclinic(3D) mode (u,v,T,S) are first bi-linearly interpolated in the horizontal direction and then a linear interpolation is used in the vertical direction from the NCOM model. The bathymetry at the open boundaries is matched between the child domain(ROMS) and the parent domain (NCOM) over a narrow transition region to avoid undesired extrapolations (Mason et al.,2010 [21] ; Pham et al., 2016 [22] ). The slowly varying signals (in the sub-tidal range) of the barotropic mode are forced by the regional NCOM model and the barotropic tide is specified using ROMS spectral forcing. The barotropic velocities,computed from the 3D velocity field, and water levels are bi-linearly interpolated in the horizontal direction and then a 21-point loess quadratic low-pass filter (40 h cutoff) is applied to filter the tide signals present in NCOM.

2.5. Open boundary conditions

Open boundary conditions are obtained from the AmSeas NCOM forecast model, with the northern and southern boundaries being located in deep water while eastern and western boundaries cross the shelf but are placed away from the shore to avoid problems due to land masking. The Chapman boundary condition (Chapman, 1985 [23] ) and the analogous Flather condition (Flather, 1976 [24] ) are used for sea surface height and barotropic momentum, respectively. In the Chapman boundary condition outgoing barotropic waves are radiated at the shallow-water wave speednd in the Flather condition deviations from exterior values are propagated out at the speed of external surface gravity waves. This combination offers a stable and robust option for the barotropic mode at the boundary, which minimize reflections of outgoing signals and on inflow conditions captures transports and water levels from the parent model. The adaptive radiation/nudging scheme of Marchesiello et al. (2005) [25] is used for the 3D velocity and tracers, with a 0.5day?1time scale for inflow condition and a 5day?1time scale for outflow determined from a sensitivity analysis.

2.6. Tide forcing

ROMS spectral forcing is used to include the tides by specifying the amplitude and phases of 11 harmonic constituents(ms4, m4, mn4, k2, s2, m2, n2, k1, p1, o1, q1), which represent diurnal (24 h), semi-diurnal (12 h) and quarter-diurnal(6 h) periods. Tidal harmonics are extracted from the latest TPXO8 global tidal solutions using the Oregon State University Tidal Prediction Software (OTPS).

2.7. Surface forcing

Surface forcing is performed by specifying directly the momentum and heat fluxes. Meteorological conditions (with the exception of wind stress) are derived from the Coastal Ocean-Atmosphere Mesoscale Prediction System (COAMPS),a 15 km resolution application of a two-way coupled oceanatmosphere model. Solar shortwave radiation, heat and freshwater fluxes are bi-linearly interpolated to the ROMS grid and specified directly at the surface. Wind stresses are computed assuming a logarithmic profile and constant sea surface roughness ofz0= 0. 003 , using winds at 10 m height taken from the National Digital Forecasting Database (NDFD), which offer quality controlled data by the National Weather Service coming from a 2 km application of the Weather Research and Forecasting (WRF) model for the region of Puerto Rico.

3.Results and discussion

3.1. Hydrographic observations

Output from the coastal ocean forecasting system is validated with remote andinsituobservations of sea surface height and ocean currents at the shelf regions. The locations of all observations points used for validation are plotted in Fig. 2 , along with color contours of the bathymetry to provide an approximate location of the shelf.

Fig. 2. Location of NOAA tide gauges (green), CariCOOS moorings with mounted Acoustic Doppler Current Profilers (blue), chosen points for High Frequency Radar time-series validation (black) and initial location of surface floating drifters for particle tracking experiment (red). Color iso-contours of depth are plotted for values of 100 m,1000 m, 2000 m, 3000 m, 4000 m and 5000 m to provide an approximate location of the the continental shelf.

The validation is performed in terms of best aggregate,such that only the first day of forecast is considered in the computed statistics. Water levels are validated at NOAA tide gauges. Sub-surface currents are validated using Acoustic Doppler Current Profilers (ADCPs), mounted on CARICOOS moorings around the adjacent coastal regions of Puerto Rico and the U.S. Virgin Islands. Surface currents are validated against surface velocity observations obtained at a 6 km and 2 km grid by the CARICOOS High-Frequency radar (HFR)network, covering most of the southwestern part of the domain.

3.2. Model data and statistics

To differentiate the ROMS code and the implementation of the operational system presented here, present results will be referred as obtained by the Puerto Rico Operational Regional Ocean Modeling System (PRO-ROMS). Similarly, AmSeas will be used to refer to the operational system application of the NCOM model used for the ROMS model initialization and boundary conditions; the forecast from AmSeas is used as benchmark for the modeled results.

The statistics used to characterize the performance of the models are the root mean square error (RMSE) and correlation coefficient (CC). The RMSE is computed as:

The correlation coefficient is defined as the covariance,or joint probability, divided by the product of the standard deviations:

Where N are the total number of sample points,fis the forecast,ris the reference or observed value and σis the standard deviation. Validation of modeled water levels with NOAA tide gauges and ocean currents with ADCP and HFR observations is performed for a period of 30 days, using the first day of forecast from January 10thto February 9th, 2017.

3.3. Water levels

Validation of Sea Surface Height (SSH) is performed at 6 NOAA tide gauge locations around the coastal waters of Puerto Rico and the U.S. Virgin Islands. A 30-day time series of the observed and modeled water levels is shown in Fig. 3 , comparing PRO-ROMS’and AmSeas’first day of forecast with the observed water levels at the tide stations plotted in Fig. 2 . We can observe that the amplitude of SSH is generally small ranging from 30 cm to 80 cm, with mixed tides on the northern coasts (e.g. San Juan) and a stronger diurnal signal on southern coasts (e.g. Yabucoa). The computed statistics show that the average RMSE and CC are approximately 5 cm and 0.9 respectively, for both PRO-ROMS and AmSeas. Validation of modeled SSH in PRO-ROMS shows good agreement with observations, with an average error below 10% and a CC of 90%. Despite its lower resolution,AmSeas predictions of SSH are similar to those obtained by PRO-ROMS, likely due to the fact that AmSeas assimilates satellite altimeter data.

3.4. Sub-surface currents

Validation of sub-surface currents is performed at middepth where the ADCP does not have significant back scattering noise from the ocean floor. A 30-day time series of observed and simulated ocean currents is shown in Figs. 4 and 5 for the eastward (u) and northward (v) components of velocity respectively. The top box on Fig. 4 shows the comparison at Ponce, located south of the island just a few hundred meters from the shelf edge ( Fig. 2 ). For this location, a RMSE of 28 cm/s and 31 cm/s is found for the eastward component of velocity for the PRO-ROMS and AmSeas models respec-tively, indicating a slight improvement in the dominant zonal current direction for this location. However, we can observe that for some time periods (e.g. days 22–27) the modeled currents by both PRO-ROMS and AmSeas models show that the current direction is opposite to the observed value. A similar trend is found from day 1 to 4 at St. John, but at this location PRO-ROMS shows significant improvement over AmSeas as shown by the RMSE of 16 cm/s and 28 cm/s for each model respectively. The large bias at certain time periods for these locations can be explained by their proximity to the shelf edge which suggests that the error is most likely caused by a poor representation of the bathymetry due to the coarser resolution of the NCOM model used by AmSeas, which at 3.6 km resolution does not allow to accurately capture the steep shelf edges in this location. Additionally, eddies observed by HFR in the southern part of the domain may not be accurately modeled due to the poor representation of the steep slopes in these areas.

Fig. 3. Best aggregate time series of modeled sea surface height by PRO-ROMS (blue), AmSeas (green) and NOAA observed (black). The NOAA locations from top to bottom are: San Juan (ID: 9755371), Mona (ID: 9759938), Magueyes (ID: 9759110), Arecibo (ID: 9757809), Yabucoa (ID: 9754228) and Charlotte Amalie (ID: 9751639) as plotted in Fig. 2 . The x-axis denotes days (January 10th to February 10th, 2017) and the y-axis denotes SSH in m.

Ocean currents at the CARICOOS Vieques buoy, exhibit strong forcing at semi-diurnal frequencies. Comparing the forecast of northward currents from PRO-ROMS and AmSeas in Fig. 5 , it is evident that AmSeas over estimates the amplitude of these signals. The phase of the currents, which show strong semi-diurnal frequencies are in very good agreement with the observations in both models. The over-estimation in the amplitude may be explained by the fact that because of the coarser resolution of the AmSeas model, less grid points are clustered in the small channels between islands. Since the forced response from the tide is approximately the same,the magnitude of the currents may be over or under specified depending on the difference in cross-sectional area of the channels between the child and parent grids. From this observation, it is expected that smaller channels are more accurately resolved by the higher resolution coastal model.

The most significant improvement is found at the San Juan buoy, where the RMSE of eastward velocity (u) shows a RMSE of 13 cm/s for PRO-ROMS while it is 28 cm/s for AmSeas, and the error is significantly reduced to less than half that of AmSeas. The close proximity of the San Juan buoy to the shore may explain the magnitude of the improvement. For this location the coastline significantly limits the northward component of velocity, while wind stress, which is mostly persistent from the East, creates strong surface forcing. The modeling study by Gonzalez [6] , found that the wind forcing significantly affects the modeled current speeds at San Juan and St. John. This reiterates the fact that wind forcing plays an important role in the current structures and ocean circulation in San Juan, and thus using winds from a 2-km WRF model results in a better local solution than the coarser resolution COAMPS model. Validations at San Juan and St. John’s ADCP’s suggest that a 2 day spin-off is long enough to significantly improve the local solution of the PRO-ROMS system over AmSeas, even without any data assimilation. Additionally, ROMS terrain-following coordinates allows the model to resolve vertical mixing in shallow shelf areas more accurately than AmSeas, and barotropic to baroclinic flow structures may differ greatly at the shelf.

Fig. 4. Best aggregate time series of modeled mid-depth eastward velocity (U) by PRO-ROMS (blue), AmSeas (green) and ADCP observed (black). The ADCP locations from top to bottom are: Ponce, Vieques, San Juan and St. John as plotted in Fig. 2 . The x-axis denotes days (January 10th to February 10th,2017) and the y-axis denotes velocity magnitude in m/s.

Fig. 5. Best aggregate time series of modeled mid-depth northward velocity (V) by PRO-ROMS (blue), AmSeas (green) and ADCP observed (black). The ADCP locations from top to bottom are: Ponce, Vieques, San Juan and St. John as plotted in Fig. 2 . The x-axis denotes days (January 10th to February 10th,2017) and the y-axis denotes velocity magnitude in m/s.

The Taylor diagram is used to summarize the overall performance for the all ADCP locations of the PRO-ROMS and AmSeas forecasts compared to the observed currents. The statistical parameters obtained for the Taylor diagrams are computed as a simple average among the 4 locations during the same time period and using the 3 h interval available for the AmSeas NCOM model in order to use the same number of sample points (N in Eqs. (1 and 2 )). Fig. 6 (a) shows the Taylor diagram for eastward currents (u) at mid-depth measured at the 4 ADCP locations. The Figure shows a RMSE of 16 cm/s and 25 cm/s for the PRO-ROMS and AmSeas forecasts respectively, yielding significant reduction of the error (about 27%) for this component of velocity. The variance of the PRO-ROMS forecast is also notably lower than that of AmSeas and almost identical to the observed currents (12 cm/s), suggesting the amplitude of the fluctuations in PRO-ROMS and the ADCP are almost identical while Am-Seas tends to over-predict these fluctuations. However, no significant improvement is achieved in the correlation coefficient between regional and coastal models.

Fig. 6 (b) shows the Taylor diagram for northward currents (v) at mid-depth measured at the 4 ADCP locations,in which we can observe similar results to the eastward currents. The figure shows a RMSE of 8 cm/s and 11 cm/s for the PRO-ROMS and AmSeas forecasts respectively, yielding significant reduction of the error (about 36%) for northward currents. The variance of the PRO-ROMS forecast is once again lower than that of AmSeas and almost identical to the observed currents (8 cm/s), which is evident from the validation at the Vieques ADCP. The main difference in the statistics of northward currents is the improvement in the correlation coefficient, which is 0.61 for PRO-ROMS and 0.54 for AmSeas. A better representation of the bathymetry and the tides result in major improvement of northward currents(v) in PRO-ROMS at Vieques and St. John ( Fig. 5 ), which exhibit significant energy at diurnal and semi-diurnal frequencies.

3.5. Surface currents

Validation of surface eastward currents (u) with the 6 km HFR at different locations is shown in Fig. 7 . The timeseries show oscillations in the tide frequency range (12–24 h period), as well as sub-tidal variability with 5–10 day periods. Comparison with AmSeas forecasts and the observed surface currents, suggests that the combination of the low frequency transports and water levels from AmSeas with ROMS spectral forcing adequately capture the large scale flows from the parent model while generating higher frequency phenomena with good agreement to HFR observations. Comparisons at the Mona Passage (HFR1 and HFR2 locations shown in Fig. 2 ) show that AmSeas overestimates the tidal amplitudes more frequently than PRO-ROMS. The RMSE for HFR1 are 19 cm/s and 25 cm/s for PRO-ROMS and AmSeas respectively. Eastward currents at HFR2 show a reduction of almost 4 cm/s in RMSE from PRO-ROMS to AmSeas, and a substantial increase in the CC from 0.16 to 0.29, suggesting a better performance by PRO-ROMS in the Mona Passage.The RMSE at HFR3 and HFR4 are 1 cm/s and 2 cm/s lower for AmSeas than for PRO-ROMS respectively, but the CC at these two locations is slightly higher for PRO-ROMS, indicating no improvement over the parent model at these locations. Note that tidal signals are more evident in the Mona Passage (HFR1 and HFR2) where PRO-ROMS shows better performance, while no significant changes are found for locations HFR3 and HFR4, suggesting that PRO-ROMS improves the forecast where tidal frequencies dominate. However, large scale flows with periods of the order of the forecast cycle (3-5days) are still constrained from AmSeas initial condition at all locations, as seen in the validations for day 15–30.

Fig. 7. Best aggregate time series of modeled surface eastward velocity (U) by PRO-ROMS (blue), AmSeas (green) and HFR observed (black). The HFR locations from top to bottom are: HFR1, HFR2, HFR3 and HFR4 as plotted in Fig. 2 . The x-axis denotes days (January 10th to February 10th, 2017) and the y-axis denotes velocity magnitude in m/s.

Validation of surface northward currents (v) with the 6 km HFR at different locations is shown in Fig. 8 . Modeled surface currents at the Mona passage show strong oscillations at the semi-diurnal frequencies, but the observed values exhibit additional phenomena with higher frequencies (less than 6 h period). At HFR1 and HFR2 the RMSE is 2–4 cm/s lower for AmSeas compared to PRO-ROMS, but both models show an overall poor performance with RMSE value of approximately 40 cm/s. This suggests that lateral (zonal) motions in the Mona Passage may be constrained by the coastlines where the tides interact with the shelf to produce more complicated signals at the surface, and the representation of the bathymetry in both models might not be accurate enough to capture these higher frequency oscillations. The steep gradients of the shelf and the smoothing performed to the original bathymetry may cause significant modifications to the ocean floor. Another plausible explanation is the wind forcing, which at 2 km resolution may still not be enough to capture the large gradients in the wake of the island located at the west of Puerto Rico due to the persistent eastern winds. Comparison at locations HFR3 yield a RMSE of 12 cm/s for PRO-ROMS and 21 cm/s for Amseas, and at HFR4 the RMSE is 7 cm/s for PRO-ROMS and 14 cm/s for AmSeas. The CC for northward surface currents at HFR3 is 0.12 for AmSeas and 0.33 for PRO-ROMS, showing significant improvement in the correlation of the forecast over AmSeas. These results show that PRO-ROMS can significantly improve the accuracy of the forecast of surface currents, especially at the shallower shelf edge location as indicated by validations at points HFR3 and HFR4. However, large scale flows associated with lower frequencies (e.g. mesoscale eddies) can significantly constraint the PRO-ROMS forecast accuracy through the initial conditions, because these flow structures take longer to evolve than the forecast cycle.

To qualitatively analyze the spatial distribution of surface currents, color contours of surface currents with superimposed velocity vectors at the Mona Passage are shown in Fig. 9 .Visual comparison with the 6 km HFR shows a good agreement of the spatial distribution of the surface currents, although the magnitude appears to be under-estimated in some regions. Validation at the near-shore with the 2 km HFR also shows good agreement in the direction and spatial distribution of the currents. Figs. 7 –9 show that PRO-ROMS can reasonably reproduce the most important features of surface ocean currents. Due to the smaller scale features in the Mona Passage, the interaction of the tides with the bathymetry and the coincidence with the meandering wake of the atmospheric circulation around the island, currents at the Mona Passage are inherently more challenging to model, as shown by the computed statistics. However, the configuration of the grid,open boundary conditions, surface and tide forcing in PRO-ROMS allow the coastal model to simulate the spatial distribution of surface currents with better accuracy than existing models.

Fig. 8. Best aggregate time series of modeled surface northward velocity (V) by PRO-ROMS (blue), AmSeas (green) and HFR observed (black). The HFR locations from top to bottom are: HFR1, HFR2, HFR3 and HFR4 as plotted in Fig. 2 . The x-axis denotes days (January 10th to February 10th, 2017) and the y-axis denotes velocity magnitude in m/s.

Fig. 9. Color contours of surface velocity magnitude with superimposed vectors of surface velocity. The top left quadrant is a snapshot of the PRO-ROMS forecast and the top right quadrant is the forecast from AmSeas for January 16, 2017 at 11:00am GMT. The bottom-left quadrant shows the observed surface velocity by the 6 km resolution HFR, and the bottom-right quadrant shows the observed surface velocity by the 2 km resolution HFR.

3.6. Particle tracking

A drifter tracking experimental campaign was organized by CARICOOS and the University of Puerto Rico to study the dispersal of larvae in the coastal regions of Puerto Rico.The experiment is useful to test particle tracking capabilities of PRO-ROMS, with the potential of developing applications for search and rescue missions, oil spill and larvae dispersal models. The experiment consisted on the deployment of 12 floating drifters with GPS mounted on them to track their position, dropped in different coastal areas around the island.The design of the floating drifters consisted of narrow tube,approximately 1m in length, which was the axis of intersection of two perpendicular planes made up of 4 sails. The locations of 4 different drifters are shown in Fig. 2 , the rest of the drifters were released at nearby areas and their trajectories are qualitatively similar.

A simple advection algorithm was used to compare the modeled trajectories by AmSeas and PRO-ROMS surface currents. To a good approximation the floating drifters or buoys are expected to behave as passive tracers and their drifting speed should be the same as the surface velocity. The algorithm takes the initial location of the drifter and uses the ocean current forecast from PRO-ROMS and AmSeas to attempt to simulate their trajectories. Ocean current velocities were interpolated in space and time to obtain the same time step for the particle simulation. In order to account for particle dispersion, 100 virtual floats were uniformly distributed over a 3 km radius and only the mean particle trajectories(i.e. average of 100 particle locations for each deployment)are considered in the results

Trajectories obtained with PRO-ROMS are closer to the drifter than those obtained with Amseas when significant oscillations are observed, as is shown in Fig. 10 (a).The initial direction of the drifter is better predicted by the PRO-ROMS forecast, and radius of the oscillations is more similar for the coastal model than the regional AmSeas model. Linear interpolation of the flow at 3 h frequency from AmSeas and overestimation of the magnitude of tidal oscillations at the Mona Passage, shown in Fig. 7 ,lead to a larger approximation of the trajectory with respect to PRO-ROMS. The higher time resolution of the PRO-ROMS forecast contribute to the better skill found in the ROMS simulations, which capture tidal oscillations more accurately than the regional AmSeas model at the coastal zones of Puerto Rico.

The average distance of the simulated Lagrangian trajectories to the actual drifter path for the 12 deployments as function of time is shown in Fig. 10 (b) for the coastal ROMS model and the regional NCOM model. The error within the first 3 h is identical, meaning that the prediction of the initial direction for both models is the same. After 15 h of deployment the error in NCOM increases almost linearly with time while the ROMS solution is approximately constant until 40 h after initial deployment. This result is a direct consequence of the resolved scales by NCOM and ROMS, and the fact that at least half of the drifters exhibit an overestimation of the oscillations observed in the NCOM trajectories, as shown for drifter D2 ( Fig. 10 (a)).

3.7. Grid sensitivity

The horizontal grid resolution is an important choice when designing a coastal ocean forecasting system, with both physical and numerical implications. In terrain-following systems such as ROMS, the sigma coordinate transformation is known to produce errors in the computation of pressure gradient terms especially in the vicinity of steep slopes and sharp density gradients (Kliem and Pietrzak, 1999 [26] ), ubiquitous in the coastal area of Puerto Rico. The use of the Shapiro filter employed in these experiments to smooth bathymetric features and reduce the pressure gradient error is known to significantly change isolated seamounts and to remove excessive volumes of steep continental shelves. This error is greatly reduced by increasing the horizontal grid resolution, which may significantly affect the energy dissipated by the shelf and consequently the response of sea surface height and currents at the shelf. From a physical standpoint, increasing the model’s horizontal resolution strengthens the vertical stratification and enhances the velocity fluctuations, especially from 3 km to 1.5 km resolution grids (Capet et al., 2008 [9–11] ). At these smaller scales, geostrophic and hydrostatic-momentum balance start to break down and submesoscale currents emerge spontaneously in coastal ocean models (McWilliams, 2017[8] )

A grid sensitivity analysis was performed using uniform horizontal resolutions of 1/30 (3.6 km), 1/90 (1.2 km) and 1/150 (0.72 km) degrees, equivalent to refinement coefficients of 1, 3 and 5 from the parent NCOM model to avoid numerical errors at the boundaries (Mason et al., 2010 [21] ; Pham et al., 2016 [22] ). Fig. 11 shows the northward barotropic current at Vieques, a small channel of approximately 10 km in width with strong tidal currents. There is a significant reduction in RMSE from the 3.6 km to the 1.2 km grid, but no significant difference is found by further increasing the resolution. The CC does not show sensitivity to the model’s horizontal resolution. Similar results are found at other ADCP locations, and therefore it is determined that the 1.2 km resolution (1/90 degree) grid offers the best compromise between model accuracy and computational efficiency.

The spectra of the measured and modeled sea surface height at the San Juan and Magueyes NOAA tide stations are shown in Fig. 12 (a) and (b) respectively. We can observe that the amplitudes at the tidal and sub-tidal frequencies are well represented at all resolutions, but significant differences are found at periods shorter than 6 h. The coarse grid with 3.6 km resolution (blue) under predicts the amplitude around the quarter diurnal (6 hr period) frequency at San Juan( Fig. 12 (a)), and increasing the resolution to 1.2 km shows that these amplitudes are better represented. This agrees well with the studies by Capet et al. (2008) [9–11] which report an important transition when increasing the resolution from 3 km to 1 km grids. The higher frequency modes at 0.8 and 2.5 cph(cycles per hour) are not reproduced by PRO-ROMS, suggesting that the bathymetry and therefore the coastline may not be well represented at these resolutions.

Fig. 10. (Left) Drifter trajectory for drifter deployment (D2). The black line represents the actual drifter trajectory, the blue line is the trajectory computed from the PRO-ROMS forecast and the green line is the computed trajectory with the AmSeas forecast. The colorbar is used to denote the time after initial deployment: the colored dots are blue at time of deployment and red 60 h after being deployed. (Right) Average distance from simulated tracks by ROMS and NCOM to floating drifters as a function of time for all 12 simulated trajectories.

Fig. 11. Best aggregate time series of modeled barotropic northward velocity at Vieques ADCP using uniform horizontal grid resolutions: 1/30 (blue), 1/90(green) and 1/150 degrees (red). The x-axis denotes days (January 10th to February 10th, 2017) and the y-axis denotes velocity magnitude in m/s.

At Magueyes a clear peak is found at the natural resonant period of Magueyes island [5] around 1.2cph, corresponding to a 50 min period. The coastal seiche is not reproduced with the 3.6 km resolution grid, but it is captured by the finer 1.2 km and 740 m resolution grids. More importantly,Fig. 12 (b) shows that there is a phase shift error associated with the model’s horizontal resolution. Both the error in the phase shift (frequency) and associated energy (amplitude) is reduced with increasing resolution. This result demonstrates the importance of horizontal grid resolution in simulating higher frequency phenomena in coastal areas where the shape of the basin affect natural resonant periods.

Fig. 12. Fast Fourier Transform of sea surface height for PRO-ROMS at 1/30 (blue), 1/90 (green) and 1/150 (red) degree resolution; (left) validation at the San Juan NOAA tide station (Station ID: 9755371) and (right) at the Magueyes NOAA tide station (Station ID: 9759110) are shown in black.

4.Conclusions

An operational Coastal Ocean Forecasting System has been developed, implemented and validated by the Caribbean Ocean Forecasting System (CARICOOS) which encompasses the region of Puerto Rico and the U.S. Virgin Islands. The Puerto Rico Operational Regional Ocean Modeling System(PRO-ROMS) is based on the ROMS ocean model, producing a 3-day forecast of ocean currents, water levels, temperature and salinity. An offline nesting approach is taken with a refinement coefficient of 3 taking initial and lateral boundary conditions from the U.S. Naval Oceanographic Office (NAVOCEANO) operational AmSeas model forecast, a 3.6-km resolution application of the regional Navy Coastal Ocean Model (NCOM) that encompasses the Gulf of Mexico and Caribbean Sea. Meteorological conditions are interpolated from the Navy’s COAMPS model with the exception of surface stresses, which are computed from a 2-km application of the WRF model used by NCEPs National Digital Forecast Database.

Validation of water levels in PRO-ROMS forecast with NOAA tide gauges show good agreement with average RMSE below 10% and CC of 0.9, but no significant improvement over the AmSeas system. Validation of sub-surface ocean currents with observed ADCP data shows good agreement with the PRO-ROMS forecast. The average root mean square error for eastward currents at mid-depth for all 4 ADCP locations is 16 cm/s for PRO-ROMS and 25 cm/s for AmSeas,and northward currents yield an average error of 8 cm/s for PRO-ROMS and 11 cm/s for AmSeas. The significant error decrease of the eastward component (u) in comparison to the northward component (v) in PRO-ROMS over AmSeas is attributed in part to the higher accuracy of the NDFD model over the COAMPS model, which uses a 2 km resolution grid compared to the 15 km resolution employed by COAMPS.Since the winds are persistent from the east, the 2 km WRF solution captures the wind curl generated by the complex topography near the coast. There is no significant improvement in the correlation coefficient for eastward transport in PROROMS compared to AmSeas, but the CC for the northward component increases from 0.54 in AmSeas to 0.61 in PROROMS. These results show an overall significant improvement in the accuracy of the modeled sub-surface currents from regional AmSeas and the coastal PRO-ROMS models, reducing the overall root-mean square error from 25 cm/s to 16 cm/s for eastward currents and a reduction from 11 cm/s to 8 cm/s for the northward currents, as well as increasing the average correlation coefficient from 0.55 to 0.61. The Taylor diagrams also show that observed currents at the ADCP locations show the same the same variance as the PRO-ROMS modeled currents, while AmSeas variance is significantly higher,especially for the eastward component.

A qualitative comparison of surface currents with the 6 km HFR show good agreement in the spatial distribution in shelf areas, as well as in the smaller scales found near-shore when comparing with the 2 km HFR observations. The time series of surface ocean currents show that PRO-ROMS adequately captures tidal frequencies with good agreement. Surface currents at the shelf edge (HFR3 and HFR4) visually show improvement in the PRO-ROMS forecast compared to AmSeas.Comparison at locations HFR3 yield an RMSE of 12 cm/s for PRO-ROMS and 21cm/s for Amseas, and at HFR4 the RMSE is 9 cm/s for PRO-ROMS and 14 cm/s for AmSeas. The CC for northward surface currents at HFR3 is 0.12 for AmSeas and 0.47 for PRO-ROMS, showing significant improvement in the correlation of the forecast with the observations.

To determine the forecast skill of simulated trajectories in PRO-ROMS and AmSeas, a particle tracking experiment was designed to test and develop applications for search and rescue missions, oil spill and larvae dispersal models. The experiment consists of 12 floating drifters with GPS mounted on them to track their position, dropped in different coastal areas around the island. The performance of the ROMS modelled tracks is improved over the NCOM tracks when significant oscillations are observed. The initial direction of the drifter is better predicted by the PRO-ROMS forecast, and radius of the oscillations is more similar for the coastal model than the regional AmSeas model. Linear interpolation of the flow at 3 h frequency from AmSeas and overestimation of the magnitude of tidal oscillations at the Mona channel. The average distance of the simulated Lagrangian trajectories to the actual drifter path for the 12 deployments as function of time show significant improvement in the forecast skill of the simulated trajectories by PRO-ROMS compared to AmSeas.

A grid sensitivity analysis was performed, using uniform horizontal grid resolutions at 1/30, 1/90 and 1/150 degrees,corresponding to 3.6 km, 1.2 km and 0.72km resolutions.Computed statistics of ocean currents and comparisons in the frequency domain via FFT of SSH show significant differences when increasing the resolution from 3.6km to 1.2km, in agreement with studies by Capet et al. 2008 [9] –[11] which report an important transition when increasing the resolution from 3 km to 1 km grids. Further increasing the horizontal resolution provide a better representation of the bathymetry and coastlines, which affect the resonant frequency at Magueyes but do reduce the model error significantly when validated with observations.

It should be noted that the improved forecast skill in the PRO-ROMS system over AmSeas is expected in shelf regions where the ADCP’s are located and where the topography plays an important role in high frequency ocean circulation processes (i.e. submesoscale currents). For example, it was shown that in small channels such as the Vieques passage and at the very near shore as is the case of San Juan, PROROMS has significantly lower root mean square error, and higher correlation coefficients. However, the assimilation of sea surface height employed by NCODA in NCOM implies that ocean processes close to geostrophic balance, such as mesoscale eddies, are expected to be better resolved by Am-Seas. Nevertheless, these evolve at timescales longer than the forecast cycle (3 days), such that the PRO-ROMS forecast is benefited from the initial conditions taken from AmSeas while resolving faster evolving coastal phenomena such as tides, which are imposed at the open boundaries.

Even though the coastal ocean forecasting system presented here does not currently employ assimilation of ocean variables, ROMS contains the tangent linear and adjoint formulations necessary for advanced assimilation techniques such as 3DVAR and 4DVAR data assimilation methods. Further improvement of the forecast may be achieved by assimilating data from ADCP and HFR measurements, and the coastal model presented here lays the foundation for such improvements.

This work was funded by CARICOOS (NOAA IOOS Grant NA11NOS0120035).