Many national weather services rely on open source forecasting models that have become available in recent decades. Due to the improvement of computational equipment and observational methods as well as the performance of numerical weather prediction models, it is now possible to provide results with good spatial resolution and increased the complexity of the simulated physics (Warner et al. 2010; B.S. Powell et al. 2009). Forecasters around the world might benefit from more information from operational near-real-time forecasting systems for their area of interest, especially if these could provide fully automatic products about the ocean-atmosphere state.
At the center of atmospheric physics at the National Institute of Meteorology of Cuba (InsMet), numerical models have been used for weather forecasting since 2006 (Mitrani et al. 2006). The MM5V3 (Fifth Generation Mesoscale Model) has been run twice per day for 72 hours of forecast with a configuretion of 3 domains (81Km, 27km, 9km) with all completely process automated and this configuretion was used in combination with the WAVEWATCH III (Third-generation numerical wave model - WW3) in order to predict waves in deep waters for cases of hurricanes (Mitrani et al. 2011). Subsequently, the SWAN (Simulating Waves Nearshore) model was included for wave prediction in shallow waters in combination with MM5V3 and WW3 (Pérez and Mitrani 2013). The modelling system presented constitutes a combination of WRF (Weather Research and Forecasting) atmospheric model, the WW3 and SWAN cited previously, as well as the ROMS (Regional Ocean Modeling System) hydrodynamic model.
This manuscript focuses on the study of Intra-Americas Sea (IAS) and Cuba which is located in the tropical area of North Atlantic. Here is presented the design of the operational forecasting system as well as statistical verification against observational data within the model domains.
WRF is a non-hydrostatic, primitive-equation, mesoscale model with advanced dynamics, physics and numerical schemes that had been widely used for many meteorological applications around the world. This model was designed to be used in operational and research activities from scales of meters to global simulations and is a result of collaboration between many institutions and university scientists (Skamarock et al., 2008). Here we use version 3.5 from which we generated the simultaneous forecast for two domains with two-way nesting allowing the interaction between the mother and child domain and centred in -80 degrees west longitude with 22 degrees north latitude. The first domain has a spatial resolution of 18 km and has 26257 grid points distributed 217 x 121 covering the Intra-Americas Sea. The second domain has 40600 grid points distributed 280 x 145 and 6 km of spatial resolution covering Cuba (see Figure 1). In both domains, the Lambert projection was employed with 30 vertical pressure levels described with sigma coordinates and was established with 50 hPa as the top level. Those domains; extensions, vertical and spatial resolutions were selected taking into account our area of interest as well as the computational power to provide results in reasonable time to forecasters. The initial and boundary global conditions were derived from the Global Forecast System (GFS) provided online by the US National Weather Service with a spatial resolution of 0.5 degrees updated every 6 hours.
The Rapid Radiative Transfer Model (RRTM) (Mlawer et al. 1997) was used for longwave radiation. It is a spectral-band and accurate scheme which use tables for efficiency. The shortwave radiation was parametrized with Dudhia scheme (Dudhia 1989) which is a simple downward integration of solar flux that uses tables for clouds. Both schemes were taken from MM5 (Skamarock et al., 2008). For the surface layer, we used the Eta similarity scheme (Janjic 1996, 2002) which is used in Eta model based on the theory of Monin and Obukhov (1954) and Zilitinkevich (1995). For the planetary boundary layer, we used the Eta operational scheme Mellor-Yamada-Janjic (Mellor and Yamada, 1982), (Janjic, 1994). The land surface parameterization used was selected according to the Noah Land Surface Model scheme operating with four layers, which improves the urban treatment and consider the surface emissivity properties (Chen and Dudhia, 2001). The microphysics scheme used was Purdue Lin, (Lin et al. 1983) which consist in a sophisticated scheme with water vapour, rain, cloud ice, cloud water, snow and graupel, appropriate for high-resolution simulations. Lastly, the Kain-Fritsch scheme is used for the cumulus parameterization (Kain, 2004) which is based on the studies of Kain and Fritsch (Kain and Fritsch 1990, Kain and Fritsch 1993), which allows deep and shallow convection.
ROMS is a free-surface, terrain-following-coordinate primitive equation ocean model that has been used for several ranges of applications by the scientific community and it is based on the Boussinesq approximation and hydrostatic vertical momentum balance (Shchepetkin and McWilliams, 2005). In this work we used the IRD (Institut de Recherche pour le Developpement) version ROMS_AGRIF (Penven P., et al 2006) (Debreu, L., et al 2011), which uses the AGRIF (Adaptive Grid Refinement in Fortran) grid refinement procedure and include a toolbox for ROMS pre- and post-processing called ROMSTOOLS (Penven P., et al 2007). This tool was employed in the creation of the computational domain containing 139264 grid points distributed at 544x256 points covering the Cuban coastal zone and extends approximately between -87.68°W and -72.4°W and between 18.26°N and 25.48°N, as shown in Figure 1. The horizontal spatial resolution was of 3 Km and 16 vertical levels using the general bathymetric chart of the oceans GEBCO (IOC, IHO and BODC, 2003) with 30 arc-second of spatial resolution.
The atmospheric forcing data was taken from the WRF in the second domain and included surface winds, long- and short-wave radiation at the surface, precipitation, surface air temperature, pressure and relative humidity using the bulk parameterization of air-sea fluxes (Fairall et al. 1996, Fairall et al. 2003). The ocean open boundaries conditions were obtained from HYCOM (HYbrid Coordinate Ocean Model), provided by the National Ocean Partnership Program and the Office of Naval Research. Data assimilative products using HYCOM are funded by the U.S. Navy. The output is publicly available at http://hycom.org in the real-time forecast with 33 vertical levels and 1/12 degree of spatial resolution. These data are horizontally and vertically interpolated on the ROMS terrain-following grid. Eight tidal constituents (M2, S2, N2, K2, K1, O1, P1, and Q1) are included in the ROMS simulations and were obtained from the Oregon State University global models of ocean tides TPXO7 (Egbert and Erofeeva, 2002). Tides were provided with 1/4 degree of resolution on a full global grid and were then interpolated to the ROMS model grid.
WAVEWATCH III (WW3) is a third generation numerical wind-wave model that has been developed at the National Center for Environmental Prediction (NCEP) (Tolman, 2002a). Version 2.22 was used here. This model uses an explicit numerical scheme, which is entirely modular, allocatable and is written in FORTRAN 90 with the possibility of being compiled for a distributed memory environment using the Message Passing Interface (MPI) (Tolman 2002b). The computational domain has a spatial resolution of approximately 6 km, and covers a region between -98.2°W and -62.8°W and between 13.1°N and 30.2°N, as is shown in Figure 1. The bathymetry for this domain was created from GEBCO using linear interpolation.
From the first domain of WRF the surface winds components data were taken to run WW3 which produce the offshore wave boundary conditions necessary as forcing for the higher-resolution nearshore wave model SWAN. In this work, an input and dissipation source term (Tolman and Chalikov 1996) was applied. The nonlinear wave-wave interactions used was the Discrete Interaction Approximation (DIA) (Hasselmann et al., 1985a, 1985b) and for bottom friction was used the empirical and linear JONSWAP (Joint North Sea Wave Project) model (Hasselmann et al., 1973).
SWAN is a third generation numerical wave model originally developed for wave simulation in shallow water and has been developed at Delft University of Technology (TU Delft) (Booij et al., 1999). This model employs implicit numerical schemes allowing increase the time steps of calculations without stability problems. In this work, we used a configuretion with 2 computational domains, each one with 1.5 km of spatial resolution, covering all the Cuban territory divided into a western and eastern zone. The western domain has 73216 grid points distributed at 416x176 points whereas that the eastern domain has 126976 grid points distributed 496x256, as is shown in Figure 1. SWAN is forced along the outer boundaries by the spectral outputs from WW3 and uses the wind fields generated within the second domain of WRF. The bathymetry was created in the same way as for the WW3, i.e. from GEBCO using linear interpolation.
In this work, we used SWAN cycle III version 40.81 with third generation and non-stationary mode, spherical coordinates and the BSBT (Backward Space Backward Time) numerical scheme. The exponential growth of wind input was used through the parameterization developed by (Komen et al. 1984). The depth-induced wave breaking was modelled using the parameterization of (Battjes and Janssen 1978) with proportionality coefficient of the rate of dissipation (α=1.0) and the ratio γ of maximum individual wave height over depth (γ=0.73). The bottom friction dissipation was modelled by the JONSWAP formulation already cited with friction coefficient cfjon = 0.067 m2 s−3, and the Komen formulation for whitecapping (Komen et al. 1984).
Figure 2 shows the workflow diagram used to develop this study. This scheme counts with global atmospheric, ocean, boundary and initial condition obtained from GFS and HYCOM respectively. Both global models are shown in Figure 2 bordered by dash lines, and their outputs are available online. In the referred diagram the shaded components are the principal components of the system.
The numerical prediction system created starts with the downloads of the initial and boundary conditions obtained from GFS. Once downloaded these files are copied to the cluster where pre-processing processes and WRF model start performing simulations. While WRF is running, the outputs are printed automatically in separate files, post-processed and graphs are displayed directly on the internal website. Once the WRF run is completed, wind fields from the first domain of WRF at 10 meters are used with WW3 to generate the wave parameters in deep waters and create the boundary conditions for the subsequent run of SWAN in the western and eastern zone of Cuba. Once the WW3 run is finished, SWAN simulations are launched in each one of the previously discussed domains using wind components obtained from the second domain of WRF. Using the HYCOM data obtained to the second domain of WRF, the ROMS model runs only once per day at 00Z. The whole mechanism is performed fully automated and currently available at InsMet.
Table 1 shows a general description of the models that constitute the numerical prediction system. The ROMS model runs only at 00Z, while the other models twice daily. All models generate forecasts for up to 72 forecast hours.
The forecast results were compared with two different sources of observations, (see the Figure 3): 31 weather stations from the national weather stations network of Cuba, provided by the climate center of the meteorological institute and 20 buoy observations obtained from the NDBC (National Data Buoy Center), which were collected and made freely available by the National Oceanic and Atmospheric Administration's (NOAA)/NDBC.
All the forecasts for each model were verified against observations (weather stations and buoys) for one year of forecast (2013) starting every day at 00z by using the bilinear interpolation from the grid-based forecasts to the observational sites. The continuous verification was done by calculating the MAE (Mean Absolute Error). This is a simple index that measures the forecast accuracy, providing the absolute value of the differences between the forecast and observation. This index does not indicate in which direction the deviations are and can take values in the range between 0 and ∞, with 0 as an optimal score. The categorical verification was made by the use of contingency tables separating the variable, in this case precipitation, in different thresholds (i.e., 0.1, 1, 2, 3, 4, 5, 10, 15, 20, 25, 30, 35, 40 mm/day). These precipitation thresholds (t) are considered to define the presence or not of a rain event and the forecasts are separated into a contingency table for each one of them. When the forecast (f) and observations (o) were ≥ t this was considered a hit. This is marked in table 2 with the letter 'a'. When f ≥ t and o < t, this was considered a false alarm. This is marked in the table with the letter 'b', when f < t and o ≥ t, these are misses and marked with the letter 'c'. When f<t and o<t, these are correct negatives, and are marked with the letter 'd'.
From the contingency table, some statistical skill measures were calculated and plotted to verify the evolution for each threshold used. The probability of detection (POD), or hit rate, is the number of hits (a) divided by the total number of events in which the observations were yes (a+c). This index gives a simple measure of the fraction of events observed that were successfully predicted by the model and can take values between 0 and 1, with 1 as the optimal score. The false alarm ratio (FAR) is the number of false alarms (b) divided by the number of events in which the forecasts were yes (a+b). This index gives a measure of how often the rain was predicted by the model but was not observed, and can take values between 0 and 1 with 0 as an optimal score.
In this paper, we also present a Taylor diagram (Taylor, 2001), which provides a summary of how the variables patterns predicted are consistent with the observations using the correlation coefficient and the standard deviation normalized for each quarter of the year.
In the following discussion, we focus our attention on the verification of the forecasting system created, using a year of observations (2013) as well as a year of numerical forecast (each day for 72 hours at 00z as the initial condition), to evaluate the performance of WRF, WW3 and ROMS. For the WRF model the temperature at 2 meters, sea level pressure, wind speed at 10 meters and precipitation, were verified. For WW3 the significant wave height was verified and in case of ROMS the sea surface temperature. The wave verification was done only for WW3 model outputs because no observational data are available for the SWAN domains.
Figure 4 shows separate quarterly and annual average evolution within 72 hours of forecasting in 2013, the mean absolute error to the temperature at 2 meters, sea level pressure, wind speed at 10 meters evaluated with weather stations in Cuba and wind speed at 10 meters evaluated with the buoys. The largest absolute error values for the temperature at 2 meters Figure 4 a) appear in the first 6 months of the year being above the annual average, while the smaller absolute errors are concentrated in the October-November-December (OND) with values below 1.5 degrees Celsius in most times of forecast. In case of the sea level pressure (Figure 4 b), the mean absolute error values for all the quarters are below 1 hPa, and this was the variable that was best forecasted by WRF in this work. The wind speed was verified against weather stations in Cuba and the buoys Figure 4 c) and d) respectively. This variable for each case of evaluation was above the annual mean values (around 2 m/s, 1.4 m/s for weather stations and buoys respectively) in the quarters OND and January-February-March (JFM), while the wind speed was mainly below the annual mean values in April-May-June (AMJ) and July-August-September (JAS).
It is essential to verify WRF wind forecasts using the buoys because it is interesting to know what the skill score introduced in wave and ocean models by the forcing from the atmospheric model is.
Figure 5 shows the evolution of POD (a) and FAR (b) for each threshold of precipitation in mm/day. The POD tends to decrease while the thresholds are increasing and the best quarter for this index is AMJ, especially for thresholds above 5 mm/day. The quarter with less POD was JFM, and this behaviour is also reflected with the FAR. Here AMJ stands out with smaller values and JFM with higher values in most of the thresholds respectively.
Figure 6 shows a separate quarterly and annual average evolution within 72 hours of forecast in 2013 of the mean absolute error of the significant wave height in the WW3 model. This variable shows the quarters OND and JFM above the annual mean values while AMJ and JAS are concentrated below the mean values. This is consistent with the wind behaviour shown already in Figure 4 c) and d), especially in this case when WW3 was forced only by the wind components from WRF. All the quarters have below 0.35 meters of absolute differences with the observations.
The ROMS model was verified by using only one buoy since this was the only one that was inside the ROMS domain and with measurements for the study period. The buoy number 42056 was selected, and the Sea Surface Temperature was verified. Figure 7 shows the quarterly and annual mean evolution within 72 hours of forecast in 2013 of the mean absolute error for this variable. The months with maximum errors are JAS with values between 0.6 and 0.7 Celsius degrees, while JFM were the best months with absolute errors between 0.25 and 0.35 Celsius degrees. An annual mean of around 0.5 Celsius degrees stands out as reasonably good for the forecast of this variable, which plays an essential role in the development and intensification of hurricanes in the study area.
Figure 8 shows the Taylor diagram for the variables temperature at 2 meters (TEMP) in red, sea level pressure (SLP) in green, wind speed in the weather stations (WSPD-S) in magenta, wind speed in the buoys (WSPD-B) in black, significant wave height (HS) in yellow and sea surface temperature (SST) in blue. Each variable is also represented for each quarter of the year. In this diagram, the standard deviations are normalized by dividing the standard deviation of the forecast variable by the standard deviation of the observations. This makes it possible to plot all the variables in the same figure. The radial coordinate represents the standard deviation normalized, and the angular coordinate represents the correlation coefficient.
The temperature at 2 meters and the significant wave height have a good correlation coefficient of around 0.9 and standard deviations slightly lower compared to the reference in all the quarters. The sea level pressure and the wind speed in the buoys have a good correlation coefficient as well, with small differences between the standard deviations for each one and the references, with a slight tendency to overestimate the reference in most of the quarters. The wind speed in the weather stations shows correlation coefficients between 0.7 and 0.9 as well as an overestimation in the standard deviation. This confirms that the wind speed evaluated from the buoys is better in absolute differences as shown in figure 4 c) and d) and also in correlation and dispersion. The sea surface temperature shows the quarters AMJ and OND with correlation coefficients around 0.9 and an underestimation in the standard deviation. The quarters JFM and JAS show low correlation coefficients between 0.3 and 0.5, with JAS standing out with an underestimation in the standard deviation values and JFM without almost dispersion.
An Ocean-Atmosphere Numerical Prediction System was developed for the Inter-American Seas and the Republic of Cuba based on the combinations of four models WRF, WW3, SWAN and ROMS, and forced with global scale models that included meteorological, hydrological and tidal data. The whole mechanism of the numerical prediction system created is performed fully automated and currently available at INSMET. The principal atmospheric and ocean variables are predicted twice daily for up to 72 hours.
In order to validate the operational models, a statistical analysis was performed for one year (2013), comparing model data with observations from weather stations and buoys. The continuous verification based on the calculation of MAE shows a reasonably good skill of the WRF model to represent the different fields, where the sea level pressure stands out with absolutes differences below 1 hPa. The temperature at 2 meters showed that the first 6 months of the year have the largest errors but all the quarters are below 1.85 degrees Celsius of differences, with annual MAE around 1.65 Celsius degree which is a good result for this variable. On the other hand, the wind speed verified with the weather stations and buoys as well as the significant wave height showed that the quarters OND and JFM lie above the annual mean values while AMJ and JAS are concentrated below the annual mean values. It is possible to see how the wind speed reduced the absolute differences in the buoys due to the non-presence of the effect of friction with the surface and how the wave model, in this case, WW3 which is mainly forced by the wind, showed the largest error in the same period. For its part, the SST predicted by ROMS showed the largest errors in JAS with values between 0.6 and 0.7 Celsius degree with an annual mean around 0.5 degrees Celsius, which is reasonably good for this variable, taking into account the important role that the SST plays in the development and intensification of hurricanes in the study area.
The accumulated precipitation in 24 hours was verified by using a contingency table for different thresholds from which the indexes POD and FAR were derived, showing AMJ as the best quarter with the largest values of the POD and smallest values of FAR for most of the thresholds selected.
We also present a Taylor diagram summarizing the general behaviour of most of the variables predicted by the system. In general, it is possible to see a good behaviour compared with the observations, with correlation coefficients above 0.7, except for the SST in the first and third quarter of the year, and with low dispersion values, except for WSPD-E for all the quarters and SST in the third quarter.
It is important to clarify that these conclusions are based on one year of simulations (2013), every day at 00Z up to 72 hours. In the future, more years will be investigated, and the statistical performance of the models will be calculated for more extended time periods. Model performance will also be compared with global models (such as GFS and ECMWF) to justify the application of high-resolution models in the forecasting system. The aim of this work was mainly to describe the selected operational setup and to show some verification aspects. At this point, it is not possible to determine whether these configuretions of domains, combinations of models, boundary conditions, and resolutions used for every component of the system are optimal, but the results are encouraging.
The operational modelling system was initially developed for the contribution to the National Meteorological Service of Cuba but will also provide data to the public, to marine activities, and to solve problems associated with coastal events, such as floods, storm surge, and oil spills.