Traditional boundaries between weather and climate prediction have been recognized as artificial (e.g., Hurrell et al., 2009; Shapiro et al., 2010) given the wide range of interactions between Earth system components on many timescales. This has motivated the seamless earth‐system prediction approach, where the same or similar earth system models are used for both short and long lead‐time forecasts (Brunet, 2010; Ruti et al., 2019). The importance of atmosphere‐ocean‐ice coupling for monthly and seasonal forecasts is well recognized (MacLachlan et al., 2015; Molteni, 2011). More recently, operational centers such as the European Centre for Medium‐range Weather Forecasts (ECMWF), U.K. Met Office, and Environment and Climate Change Canada (ECCC) are also finding benefits in ocean coupling for lead times as short as a few days (Smith, 2018), which further supports the seamless approach to environmental prediction.
Along these lines, the U.S. Naval Research Laboratory (NRL) has developed a new global coupled data assimilation and prediction system, the Navy Earth System Prediction Capability (Navy‐ESPC), with the goals of both improving short‐term forecasts and extending forecasts to multiweek timescales. The Navy‐ESPC effort is part of the U.S. National ESPC effort, which is a five‐agency effort to coordinate and accelerate the national numerical environmental prediction capability (Carman, 2017). As noted in that article and others, there are many needs for environmental predictions over various timescales. U.S. Navy applications that would benefit from skillful environmental predictions include the planning and conducting of Navy ship sorties on timescales less than a week, ship routing and prepositioning on timescales out to a month, and humanitarian assistance planning and force deployment management on monthly to seasonal timescales. Unlike many other forecasting centers, the U.S. Navy is interested in high‐fidelity forecasts of the ocean itself, including the ocean interior. As described below, this need for accurate ocean forecasts drives Navy‐ESPC design, and the high‐horizontal resolution ocean and sea‐ice components of Navy‐ESPC distinguishes it from other global coupled prediction systems.
Navy‐ESPC is built upon existing atmosphere and ocean/sea ice operational systems that are being used at Fleet Numerical Meteorology and Oceanography Center (FNMOC). The atmospheric model is the NAVy Global Environmental Model (NAVGEM) (Hogan, 2014), and the ocean and sea ice models are part of the Global Ocean Forecasting System version 3.1 (GOFS 3.1) (Metzger et al., 2014). Navy‐ESPC is unique compared to many other global coupled modeling systems (Infanti & Kirtman, 2016; Lin et al., 2016; Saha, 2014; Sun et al., 2018; Vitart et al., 2017; Zhu, 2018) in that the ocean model horizontal resolution is eddy resolving in both the ensemble (i.e., probabilistic) and deterministic configurations. In many operational coupled modeling systems, the ocean and sea ice components are included primarily to provide better atmospheric bottom/surface boundary conditions but only at eddy‐permitting resolution (~1/4° or coarser). However, users of Navy‐ESPC products need high‐resolutions predictions from each component.
A high‐resolution ocean model is needed to simulate ocean energetics (Chassignet & Xu, 2017; Hogan & Hurlburt, 2000; Thoppil et al., 2011) as the Rossby radius of deformation is much smaller in the ocean than in the atmosphere (Feliks, 1985; Hallberg, 2013). Thoppil et al. (2011) compared ocean modeling results from a 1/12° and 1/25° ocean and found that the Eddy Kinetic Energy (EKE) in the 1/25° nonassimilative simulation and 1/12° data assimilative reanalysis matched observations more closely than the 1/12° nonassimilative simulation. Chassignet and Xu (2017) examined Gulf Stream separation in ocean models ranging from 1/4° to 1/50° and determined that the Gulf Stream was simulated most accurately in the 1/50° ocean. In examining the impact of high‐resolution ocean modeling on atmosphere predictions, Hewitt (2017) determined that the high‐resolution ocean was beneficial to the atmosphere near western boundary currents.
The main objective of this paper is to describe the various components of the Navy‐ESPC system and highlight key metrics used in the validation and verification process as it was transitioned to FNMOC. Due to resource limitations, the changes between the coupled system and the current operational systems could not be tested independently. The purpose of this article is to provide a detailed description of the new system, which will serve as a reference for future hypotheses‐driven papers and other studies that may use the system output. The modeling components, coupling technologies, and data assimilation in sections 2 and 3; forecasts used for the diagnostics in section 4; results in section 5; and conclusions and future directions in section 6.
The global atmosphere component of Navy‐ESPC is NAVGEM (Hogan, 2014) and the ocean and sea ice components are the HYbrid Coordinate Ocean Model (HYCOM) (Bleck, 2002) and Community Ice CodE (CICE) (Hunke &> Lipscomb, 2008), respectively (Figure 1). HYCOM coupled to CICE currently runs at FNMOC as stand‐alone GOFS 3.1. NAVGEM and HYCOM are developed at NRL and CICE is from the Los Alamos National Laboratory. These models, which are also referred to as components, are coupled together using Earth System Module Framework (ESMF) tools, in conjunction with the National Unified Operational Prediction Capability (NUOPC) layer (Theurich, 2016).
Enlarge this image.1 Figure. Schematic of the Navy‐ESPC modeling system. The atmospheric, ocean, and sea ice model components are the NAVy Global Environmental Model (NAVGEM), the HYbrid Coordinate Ocean Model (HYCOM), and the Community Ice CodE (CICE), respectively. HYCOM and CICE are part of the operational Global Ocean Forecasting System version 3.1 (GOFS 3.1). To couple the models, the Earth System Module Framework (ESMF) tools are used in conjunction with the National Unified Operational Prediction Capability (NJUOPC) framework. The colored arrows represent variables that export and import from each system. Note the land model is imbedded inside of NAVGEM.
NAVGEM (Hogan, 2014) is the Navy's global weather prediction system, run operationally at FNMOC and is used as the atmospheric component in Navy‐ESPC. The forecast model uses both grid point and spectral (i.e., spherical harmonic) representations to perform the forecast. Grid point calculations are performed for the Semi‐Lagrangian (SL) advection and in all physical parameterizations. Calculations in spectral space are performed for the Semi‐Implicit (SI) corrections to the divergent component of the winds, virtual potential temperature, and surface pressure. The dynamical core of NAVGEM is a three‐time level, SL/SI numerical integration of the hydrostatic equations of motion and the first law of thermodynamics. The SL/SI formulation was first implemented at ECMWF based on the work of Ritchie (1987, 1988, 1991) and Ritchie et al. (1995). The major difference between the Ritchie and NAVGEM formulations is that NAVGEM uses potential temperature while ECMWF uses temperature.
The Navy‐ESPC uses two different resolutions: for NAVGEM, the high resolution deterministic (Navy‐ESPCDET) is T681L60 (~18 km horizontal resolution) and the ensemble resolution (Navy‐ESPCENS) is T359L60 (~37 km horizontal resolution) (Table 1). The dynamics framework utilized for the Navy‐ESPCDET is slightly different from the Navy‐ESPCENS configuration. NAVGEM in Navy‐ESPCENS followed then operational stand‐alone NAVGEM 1.2 (T359L50) in retaining the virtual potential temperature as the prognostic temperature variable, and using the full Gaussian grid. The adiabatic correction option introduced in NAVGEM 1.2 is turned off for the ensemble resolution. The adiabatic correction option in NAVGEM is a method to address the lack of conservation of mass‐weighted virtual potential temperature, which occurs largely in the stratosphere in NAVGEM. Mass loss/gain is added/subtracted back into model. NAVGEM in Navy‐ESPCDET followed then operational stand‐alone NAVGEM 1.4 (T425L60) and uses the three‐time level SL/SI dynamical core, perturbation virtual potential temperature for the prognostic temperature variable, and uses the reduced Gaussian grid. A 60‐level vertical grid was used in both horizontal resolutions in Navy‐ESPC. This is the same vertical resolution used in the recent operational version of NAVGEM.
TableConfigurations of the Components/Models in the High Resolution Deterministic (Navy‐ESPCDET) and Ensemble (Navy‐ESPCENS)| Navy‐ESPC configuration | Navy‐ESPCDET | Navy‐ESPCENS | |
| NAVGEM | Resolution | T681 60 Levels | T359 60 Levels |
| Parameterization Differences |
Base Version: NAVGEM 1.4 (resolution T425L60) Prognostic Temperature Variable: perturbation virtual potential temperature Grid Type: reduced Gaussian grid Adiabatic Correction: Yes |
Base Version: NAVGEM 1.2 (resolution T359L50) Prognostic Temperature Variable: virtual potential temperature Grid Type: Gaussian Adiabatic Correction: None |
|
| HYCOM | Resolution | 1/25° tripole 41 layers | 1/12° tripole 41 layers |
| Parameterization Differences | Tides: M2,S2,K1,O1,N2 | Tides: none | |
| CICE | Resolutions |
1/25° tripole 4 ice categories |
1/12° tripole grid 4 ice categories |
| Parameterization Differences | No differences | No differences |
Note. NAVGEM versions in Navy‐ESPCDET and Navy‐ESPCENS both use the physics changes listed in.
For the implementation of NAVGEM into NAVY‐ESPC, we developed a new suite of model physics referred to as Coupled Version Physics (CVP) (Table 2). The CVP physics have been tailored for the coupled system with two primary objectives: (1) to improve the representation of the Madden‐Julian Oscillation (MJO), which is recognized as a key contributor to extended range predictability (Waliser, 2009), and (2) to improve the consistency between surface fluxes computed in NAVGEM and in HYCOM. The MJO is a mode of atmospheric variability in the tropics with a timescale of 20–100 days and is one of the few phenomena affecting global predictability on extended‐range timescales (Lim et al., 2018). The phase of the MJO is known to have substantial impact on many other phenomena such as tropical cyclone genesis and frequency (Camargo et al., 2008; Kim et al., 2008; Sobel & Maloney, 2000), atmospheric rivers (Bond & Vecchi, 2003; Guan et al., 2012; Mundhenk et al., 2016), and even the atmospheric circulation in the Arctic and Antarctic (Flatau & Kim, 2013; Henderson et al., 2014, 2018; Lin et al., 2009). As such, it is important for extended‐range forecasting systems to be capable of realistic simulations and accurate predictions of the MJO.
TableSummary of NAVGEM Physics Changes Between NAVGEM v1.4 and the NAVGEM Used in the Navy‐ESPC system| Parameterization/scheme | NAVGEM v1.4 | NAVGEM in Navy‐ESPC |
| Convection parameterization | SAS (Moorthi et al., 2001) | Modified Kain‐Fritsch |
| Boundary layer scheme | Louis et al. (1982) | COARE (Kara et al., 2005) |
Note. The main changes between the two atmospheric models include changes in the convection parameterization and boundary layer scheme. The Modified Kain‐Fritsch is based on Ridout et al. (2005).
To improve the representation of the MJO in Navy‐ESPC, we replaced the operational NAVGEM Simplified Arakawa‐Schubert (SAS) convection scheme (Han & Pan, 2011) with a modified version of the Kain‐Fritsch scheme (Kain & Fritsch, 1990, 1993). The new scheme is an extension of Ridout et al. (2005) and incorporates both turbulence and dynamically forced modes. As in the 2005 version, the scheme includes a modified closure formulation for the Kain‐Fritsch dynamically forced mode based on an assumed quasi‐balance of updraft parcel buoyancy at the cloud base level. This closure formulation requires the scheme to be called at every time step to adjust the cloud base mass flux in a similar manner as in the Emanuel convection scheme (Emanuel, 1991; Emanuel & Zivkovic‐Rothman, 1999). For the turbulence‐forced mode of convection, triggering by boundary layer plumes is represented in part using the mixed‐layer Richardson number convective trigger formulation described by Ridout and Reynolds (1998). Other added enhancements for Navy‐ESPC include a treatment of convective momentum transport and a modified representation of the rate of updraft‐environment mixing adapted from Peng et al. (2004).
In addition to the convection scheme changes, the other update for the CVP physics is in making the fluxes more consistent between the atmosphere and ocean. The HYCOM surface flux scheme is the Kara et al. (2005) adaptation of the COARE 3.0 scheme of Fairall et al. (2003) and is incorporated in NAVGEM within Navy‐ESPC. The coupled system framework enables for the first time in NAVGEM the inclusion of surface currents in the surface flux computation. Navy‐ESPC versions of NAVGEM and HYCOM also account for depression of the saturated vapor pressure due to salinity using the correction of Sud and Walker (1997). See Figure 2 for the list of variables coupled to NAVGEM.
Enlarge this image.2 Figure. List of fields passed to and from the model components (i.e., variables that are coupled). The color represents the component/model the field is from and the arrows represent if the fields are going into or out of the component. The long names for the atmospheric fields include: U @ 10 m = zonal winds at 10 m; V @ 10 m = meridional winds at 10 m; Temp @ 2 m = air temperature at 2 m; Spec Hum @ 2 m = specific humidity at 2 m; Density @ 2 m = air density at 2 m; Temp @ Sur = temperature of the surface (ground, ice, ocean); SW down = downward surface shortwave radiative flux; LW down = downward surface longwave radiative flux; Precipitation = total precipitation; MSLP = mean sea level pressure. The long names for the ocean variables include the following: SST = sea surface temperature; SSU = zonal sea surface current; SSV = meridional sea surface current; SSS = sea surface salinity. The long names for the sea ice variables include the following: SIC = sea ice concentration; X Ice Stress @ B = x direction ice stress at sea ice basal; Y Ice Stress @ B = y direction ice stress at sea ice basal; SW Flux @ B = shortwave radiative flux at sea ice basal; LW Flux @ B = longwave radiative flux at sea ice basal; Salt Flux @ B = flux of salt at sea ice basal; Water Flux @ B = flux of fresh water at sea ice basal; Ice U = zonal velocity of sea ice at the surface; Ice V = meridional velocity of sea ice at the surface; Ice Thick = sea ice thickness; Ice Albedo = sea ice albedo; Ice Temp = temperature of sea ice at surface.
The CVP has contributed to an improved representation of the MJO with respect to stand‐alone NAVGEM. Stand‐alone NAVGEM's representation of the MJO was not competitive compared to other global models, as shown in Jiang (2015). With the implementation of the CVP and coupling, the Navy‐ESPC model compares favorably in this regard to other state‐of‐the‐art forecast as seen through evaluation of the 17+ years of forecasts in Janiga et al. (2018) and Kim et al. (2019). These 17+ years reforecasts were part of the National Oceanic and Atmospheric Administration (NOAA) Subseasonal eXperiment (SubX) (Pegion, 2019). The CVP physics contributes to this notable improvement, as does the atmosphere‐ocean coupling in Navy‐ESPC. Figure 3 displays Hovmöller diagrams of precipitation during the Dynamics of the Madden‐Julian Oscillation (DYNAMO) observation campaign (Johnson & Ciesielski, 2013) for operational stand‐alone NAVGEM v1.4, NAVGEM v1.4 with the CVP option, and Navy‐ESPC. As shown, operational NAVGEM v1.4 did not represent the first or second MJO for this time period. With addition of the CVP physics, precipitation for this period is more similar to the observations than without the CVP physics, and including coupling within Navy‐ESPC further improves the representation of the MJO precipitation.
Enlarge this image.3 Figure. Hovmöller diagram of precipitation from NAVGEM v1.4 with NAVGEM physics, NAVGEM v1.4 with the coupled physics, Navy‐ESPC, and The Tropical Rainfall Measuring Mission (TRMM) from 1 November 2011 to 10 December 2011. Precipitation is averaged from 5°S to 5°N and the Hovmöller diagram shows precipitation from 40°E to 140°E. NAVGEM and Navy‐ESPC output represent long nonassimilative forecasts starting on 1 November 2011. NAVY‐ESPC is configured as a deterministic Navy‐ESPCENS run with 50 vertical levels.
The HYCOM is the ocean model in Navy‐ESPC. HYCOM is a primitive equation ocean general circulation model capable of nowcasting and forecasting the three‐dimensional temperature, salinity, and current structure of the global ocean. It employs potential density referenced to 2,000 m and includes the effects of thermobaricity (Chassignet et al., 2003). HYCOM's grid is uniform cylindrical from 78.64°S to 66.0°S, on a Mercator projection from 66.0°S to 47°N and curvilinear north of this as it employs an Arctic dipole patch where the poles are shifted over land to avoid a singularity at the North Pole. This grid is referred to as a tripole grid.
HYCOM uses a mixture of vertical coordinates, hence the HYbrid as part of its name. Vertical coordinates can be (1) isopycnals (density tracking), often the best coordinate in the deep stratified ocean; (2) levels of equal pressure (nearly fixed depths), best used in the mixed layer and unstratified ocean; and (3) sigma‐levels (terrain following), often the best choice in shallow water. HYCOM combines all three approaches by choosing the optimal distribution at every time step. The model makes a dynamically smooth transition between coordinate types by using the layered continuity equation. The hybrid coordinate extends the geographic range of applicability of traditional isopycnic coordinate circulation models toward shallow coastal seas and unstratified parts of the world ocean. It maintains the significant advantages of an isopycnal model in stratified regions while allowing more vertical resolution near the surface and in shallow coastal areas, hence providing a better representation of the upper ocean physics.
HYCOM is configured with options for a variety of mixed layer submodels (Halliwell, 2004) and here uses the K‐Profile Parameterization (Large et al., 1994). A more complete description of HYCOM physics can be found in Bleck (2002). The equatorial horizontal resolution of HYCOM is 1/25° in Navy‐ESPCDET and 1/12° in Navy‐ESPCENS, and both configurations have 41 layers (Table 1). The number of vertical layers is a compromise between accuracy and computational efficiency when run in an operational setting.
HYCOM versions in Navy‐ESPCDET and Navy‐ESPCENS are identical except Navy‐ESPCDET turns on astronomical tidal forcing that generates internal tides (Arbic et al., 2018) and the Navy‐ESPCENS does not. HYCOM implementation in stand‐alone operational GOFS 3.1 is very similar to that in Navy‐ESPC, and the main difference is that the Navy‐ESPC HYCOM uses the coupled sea level pressure forcing instead of constant sea level pressure forcing as used in GOFS 3.1. This change in sea level pressure forcing is needed for the internal tides in Navy‐ESPCDET.
Version 4.0 of the Los Alamos‐developed CICE model (Hunke & Lipscomb, 2008) is the sea ice component of Navy‐ESPC. CICEv4 is also currently used in stand‐alone operational GOFS 3.1. CICE includes sophisticated ice thermodynamics such as multiple ice thickness layers, multiple snow layers, and the capability to forecast multicategories of ice thickness according to World Meteorological Organization definitions. In addition, CICE has several interacting components including a thermodynamic model that computes local growth rates of snow and ice due to snowfall, vertical conductive, radiative, and turbulent fluxes; a model of ice dynamics that predicts the velocity field of the ice pack based on a model of the material strength of the ice; a transport model that describes advection of the areal concentration, ice volumes, and other state variables; and a ridging parameterization that transfers ice among thickness categories based on energetic balances and rates of strains. The CICE horizontal grid for Navy‐ESPCDET and Navy‐ESPCENS is the same as HYCOM (Table 1) for conservation of variables when coupling. There are no differences in CICE physics between the deterministic and ensemble configurations of Navy‐ESPC. A more recent version of CICE (i.e., version 5 or 6) was not used in the Navy‐ESPC code simply because the code needed to be locked for testing to meet timelines for transition to operations.
Navy‐ESPC uses the NUOPC tools wrapped around ESMF for coupling. This framework allows the components to dynamically interact at a specified coupling time step. A 1‐hr coupling time step is implemented in Navy‐ESPC. We tested a more frequent time step; however, the results were not improved enough to warrant the additional computational cost. As the analyzed ocean and sea‐ice fields are held constant during the stand‐alone NAVGEM integrations, the coupled system holds the promise of improved surface forcing for the atmosphere, particularly for forecast lead times greater than a few days. Figure 2 displays the coupling variables in Navy‐ESPC. Each model component performs its own flux calculations. The implementation of a single coordinated flux calculation was considered, but found to be impractical for this initial rollout of Navy‐ESPC with time constraints set by our operational partners. NAVGEM's use of a fully implicit mixing scheme increases the computational demands in this context, as does the large mismatch between the high‐resolution ocean grid and the lower‐resolution atmosphere grid, which results in an inefficient interpolation of variables. Mixing in NAVGEM's scheme is carried out through the depth of the atmospheric column in conjunction with the computation of surface fluxes. Surface exchange coefficients are computed based on Kara et al. (2005), and then a tridiagonal system is solved to obtain the mixing throughout the column, including the surface fluxes. To perform this calculation in the mediator, three‐dimensional fields would need to be interpolated onto the high‐resolution mediator grid, which is the same as the HYCOM grid, and ideally be performed at every NAVGEM time step (5–6 min). With the mismatch between the high‐resolution HYCOM grid and the NAVGEM grid, this interpolation was deemed to be too computationally expensive and memory intensive. The lack of a consistent flux calculation between components results in the loss of conservation. However, note that our forecasts do not go beyond 45 days, and our experience with stand‐alone NAVGEM has shown that the implicit approach provides a skill advantage for numerical weather prediction, one that we wish to retain.
Another difficulty in the interpolations is that the ESMF remapping routines initially led to grid points with undefined values near coastlines. In Navy‐ESPC, the land/ocean interpolation boundary is complicated thanks to differences in coastline representation in NAVGEM and HYCOM driven by (1) differences in horizontal resolution and (2) the definition of the land or ocean point. For example, the HYCOM land/ocean boundary is defined by a depth in the ocean in which HYCOM will begin calculations, and the NAVGEM land/ocean boundary is defined where surface fluxes are greatly altered by the surface state, which is generally at a shallower ocean depth than HYCOM's land/ocean boundary. With the ESMF remapping routines, unassigned points occurred in the HYCOM to NAVGEM and NAVGEM to HYCOM interpolations, and we developed a routine to address this shortcoming. First, we find the grid points that do not have values (from the default ESMF scheme), and call these “hole cells.” Next, we create a stencil of grid points around the hole cell and assign each stencil point as land or ocean. If the hole cell is ocean, the ocean points in the stencil are used to interpolate a value to the hole grid cell. If the hole cell is land, land points in the stencil are used for the interpolation. In our new routine, the stencil size and the number of times to search for hole cells are user defined and is only used at the start of the model simulation. A simple implementation of this algorithm has been added into version 8 of ESMF called ESMF_EXTRAPMETHOD_CREEP.
The Navy‐ESPC system utilizes two mature assimilation code bases that have been developed for GOFS 3.1 and NAVGEM forecast systems separately: the Navy Coupled Ocean Data Assimilation (NCODA) for the ocean and sea‐ice, and the Naval Research Laboratory (NRL) Atmospheric Variational Data Assimilation System Accelerated Representer (NAVDAS‐AR) for the atmosphere.
NCODA is a three‐dimensional, multivariate (3DVAR) Data Assimilation (DA) scheme (Cummings, 2005; Cummings & Smedstad, 2013) for the ocean/ice variables temperature, salinity, geopotential, vector velocity components, and sea ice concentration. All observations are processed at FNMOC through a quality control system that checks for gross error and assigns each observation a likelihood of error based on comparison to climatology and ocean analyses. Observation errors are assigned by observation type and include both an instrument error and a representativeness error that account for the observation spatial sampling and the model grid resolution. Sea surface height anomaly (SSHA), sea surface temperature (SST), and ice concentration observations are pooled into superobservations at the resolution of the model grid before assimilation. Along‐track SSHA observations from available altimeters referenced to a 1993–2015 mean height are used to derive synthetic subsurface temperature and salinity profiles via the Improved Synthetic Ocean Profile (ISOP) methodology of Helber et al. (2013), and these are assimilated together with satellite and in situ SSTs from ships, and moored and drifting buoys. Profile temperature and salinity observations come from Argo profiling floats, XBTs, CTDs, gliders, and marine mammals. Satellite sea ice concentration is assimilated in the ice analysis. Ocean data assimilation windows are typically longer than in atmospheric application, and subsurface ocean observations provide useful information for days after the observation time. The data windows for profile and SSHA observations may be extended back in time both to capture late‐arriving observations and to include a larger set of observations more capable of capturing the ocean mesoscale structure. While the atmosphere is “data rich,” the ocean is “data poor” with regard to having enough subsurface observations to accurately constrain the smaller scale (than the atmosphere) mesoscale features. We use a 5‐day window for in situ profile observations and a 4‐day window of satellite SSHA observations. For satellite and in situ SST and satellite ice concentration, the data window extends back only to the most recent observation assimilated previously. We do not presently apply bias corrections to the satellite SSTs. SST bias corrections may be explored in the future as these corrections have been shown to aid in surface flux estimates (While & Martin, 2019). The data are compared against a time dependent background field using the First Guess at Appropriate Time (FGAT) method (Cummings & Smedstad, 2014; Fox et al., 2002).
NCODA can be run in stand‐alone mode but here is cycled with HYCOM and CICE to provide updated initial conditions for the next model forecast using an incremental analysis update (IAU) procedure (Bloom et al., 1996). In an uncoupled model (i.e., GOFS 3.1), the first guess is generated by the ocean model forecast from the end of the last DA window to the end of the current DA window.
NAVDAS‐AR (Rosmond & Xu, 2006; Xu et al., 2005) has the ability for both weak and strong constraint variational assimilation, and is formulated in the terms of dual variables (i.e., observation space), the dimension of which is generally much smaller than the corresponding state (i.e., model) space. This restriction of the observation space is accomplished by discarding the unobservable degrees of freedom in the system (Bennett, 2002).
NAVDAS‐AR can process well over 100 million observations in every 6‐hr data assimilation window. After quality control and data thinning it assimilates approximately 4 million observations to create the final analysis. Within the 4DVAR data assimilation cycle, NAVDAS‐AR approximates the full nonlinear atmosphere in the Navy‐ESPC coupled forecast with a simplified version of the stand‐alone NAVGEM system to allow for a computationally efficient solution of the assimilation problem.
Navy‐ESPC implements a weakly coupled formulation of the coupled data assimilation system with two separate windows for the atmosphere and ocean. In a weakly coupled data assimilation, the coupled model forecast is used as a first guess for the independent data assimilation systems. Given this coupled first guess, the ocean (NCODA) and atmospheric (NAVDAS‐AR) data assimilation systems compute the assimilation increment independently, ignoring any potential cross‐correlations between ocean and the atmosphere. To simultaneously accommodate the 6‐hr DA cycle that NAVDAS‐AR has been designed for and the 24‐hr DA cycle that has proven optimal NCODA, a new dual cycle DA has been developed (Figure 4, Table 3).
Enlarge this image.4 Figure. Schematic of the weakly coupled DA system for Navy‐ESPC. The red arrows represented the Navy‐ESPC first guess/forecast from the coupled model. The yellow bars represent the observational window used for NAVDAS‐AR and the blue bar represents the ocean and sea ice observational window used for NCODA. The light red boxes represent the 3‐hr IAU period when the increments produced by NCODA and NAVDAS‐AR are inserted into the Navy‐ESPC coupled forecast. The yellow and blue arrows pointing down from the coupled forecasts represent the times of the forecast that are used for first guess for (yellow) NAVDAS‐AR and (blue) NCODA. Note that NCODA appends the forecast together at 00Z and 06Z for its first guess. Also note that sea ice increments are directly inserted at 09Z.
| NAVDAS‐AR | NCODA | |
| Stand‐Alone and Weakly Coupled Implementation | ||
| Forecast Model(s): | NAVGEM | HYCOM and CICE |
| Method: | 4DVAR Hybrid | 3DVAR FGAT |
| Update Window: | 6 hr (±3 hr) | 24 hr (±12 hr) |
| Stand‐Alone Implementation | ||
| Updating Method: | Direct Insertion | 3‐hr IAU |
| Weakly Coupled Implementation | ||
| Updating Method: | 3‐hr IAU | 3‐hr IAU |
The centers for the 6‐hr atmospheric DA widows are on 00Z, 06Z, 12Z, and 18Z. Coupled forecast from the previous update cycle is performed to at least 3 hr after the center point of the next update cycle, with the 6‐hr section of the forecast centered on the next analysis time used as the first guess to initialize the atmospheric portion of the coupled model for the subsequent center point of the next update cycle. For example, from the 00Z coupled forecast, we use forecast hours valid between 03Z and 09Z as the first guess for the initialization at 06Z. Next, we run NAVDAS‐AR on this first guess field to obtain the analysis increments/corrections valid at 06Z. The correction is then inserted into the coupled model by IAU over the 3 hr prior to the analysis time. For example, for the 06Z centered window, the coupled forecast starts at 03Z with the increment incorporated over the 03Z to 06Z using IAU. Thus, 12‐hr coupled forecasts are needed for each atmospheric DA cycle, 3 hr needed for IAU plus the 9 hr needed for the next first guess. The 3‐hr IAU in NAVDAS‐AR is a new capability implemented specifically for the Navy‐ESPC as the FNMOC operational NAVDAS‐AR implementation does not use an IAU. The forecasts started at 06Z are run 12 hr longer than the forecasts started at 00Z, 12Z, and 18Z to provide a first guess for the 24‐hr ocean DA window centered at 12Z, described in the next paragraph. Unlike ECMWF's system (Laloyaux et al., 2016), but similar to our stand‐alone atmospheric system, Navy‐ESPC uses only a single outer‐loop in the atmospheric NAVDAS‐AR system.
The NCODA ocean DA window center is at 12Z and uses a 24‐hr window that is 12 hr before and after the 12Z center (i.e., 0 to 24Z) with an extended look backward for late arriving observations as described above. The first guess for NCODA combines the 0 to 6 hr coupled forecast center at 00Z, and the 0 to 18 hr coupled forecast at 06Z, which corresponds to times 06Z to 24Z. As in the atmosphere, a 3‐hr IAU is used to update the ocean at the 12Z center point, which corresponds to the coupled model forecast starting 3 hr before 12Z, at 9Z. For the sea ice update, the analysis produced by NCODA is directly inserted into CICE at 09Z without any incremental insertion. The ocean does not update with increments/corrections at 00Z, 06Z, and 18Z.
Using the IAU method for both the atmosphere and ocean leads to the coupled model transferring some of the initialization information from the atmospheric observations to the ocean and vice versa during the IAU period. Since both the ocean and atmosphere are updated at 12Z, only the forecasts initialized at 12Z are used for the forecast verification.
To manage the complexity of the weakly coupled data assimilation system, such as task dependencies, we rely on the Cylc workflow manager (Oliver, 2019). The workflow developed for Navy‐ESPC represents the combination of workflows developed for each standalone component, with the coupled forecast model as a common component. Additionally, the ensemble configuration makes use of independent workflow execution for each ensemble member.
To generate initial conditions for the subseasonal ensemble forecast system, we employ the ensemble of data assimilations approach (EDA) (Houtekamer et al., 1996). In our implementation, one member is used as a control member and is cycled using the ensemble‐resolution version of the Navy‐ESPC system. In addition to the control member, we also cycle 15 perturbed members, which assimilate “perturbed” observations that included Gaussian noise errors consistent with the observational error statistics prescribed by the data assimilation system. The number of ensemble members was determined by the computer resources available at FNMOC. In NAVDAS‐AR, these random errors are added to the innovation vector after thinning and gross quality control checks but before the bias correction is applied to the radiances. The NAVDAR‐AR system uses a Variational Bias Correction (VarBC), method which adds the bias coefficient predictors to the control vector for the data assimilation, solving for them at each cycle point (Auligne et al., 2007; Dee, 2005). For NCODA, the raw observations are perturbed before NCODA generated synthetic profiles based on sea surface height anomalies and before thinning and superobservations are computed. Due to the relative numbers of real and synthetic subsurface observations, the largest impact of the ocean observation perturbations is through the ISOP profiles generated by perturbing the SSHA and SST values use as inputs to the ISOP system. This generates observation spreads consistent with the measured model forecast errors (not shown). Although the perturbation of in situ profile observations has limited impact in areas regularly sampled by the altimeters, a method is being developed that adjusts the profile vertically to fit a random perturbation of the profile steric height, as is the simpler method of perturbing the profile location on the model grid after the profile selection process is complete. Either of these should improve on the present approach of perturbing profile level observations independently with no accounting for correlated error.
Unlike traditional EDA systems (El Ouaraini & Berre, 2011; Houtekamer et al., 1996), the Navy‐ESPC system does not feed the ensemble variance back into the covariance model used for data assimilation. Furthermore, we do not center the ensemble system on either the control member or the high‐resolution analysis from the Navy‐ESPCENS system. We made these choices because of the extreme computational expense of the system that made it impractical to synchronize all the ensemble members at every update cycle.
Prior research (Bowler, 2017) indicated that EDA systems tend to be underdispersive at initial time with sluggish growth of unstable baroclinic modes in the atmosphere. While fully realizing this drawback of the pure EDA system, we chose not to include either stochastic noise or dynamically condition posterior ensembles through inflation or addition of atmospheric singular vectors. This choice was motivated by a compressed development timeline for the first version of the system, and we plan to address these shortcomings in upcoming versions through methods for representing the stochastic effect of unresolved processes (Crawford et al., 2020).
For this paper, we run Navy‐ESPC with two different configurations/resolutions (Table 1). For each configuration, reanalysis cycles started on 15th December 2016 at 12Z. The atmospheric Initial Conditions (ICs) for the starts of the reanalyses were taken from the operational NAVGEM for both configurations. For HYCOM and CICE, the Navy‐ESPCENS ICs are from operational stand‐alone GOFS 3.1 and the ICs for Navy‐ESPCDET are from a 1/25° stand‐alone GOFS 3.5 (Metzger et al., 2019) that has been approved for operational testing at FNMOC.
For the ensemble configuration, the ensemble reanalysis started from the same initial condition on 15th December 2016. We allowed our EDA system to develop spread in the ensemble initial conditions by cycling the system for 1.5 months (until 1st February 2017). Preliminary results suggested that the ensemble spread in the EDA system stabilized after 1–3 days in the atmosphere and after 1–2 month in the ocean (for most variables). Some small increase in the sea surface height spread continued in the Western boundary regions up to 6 months into the reanalysis.
The reforecasts results presented in this paper for a single member of the Navy‐ESPCDET are from the period 1 January 2017 to 31 December 2017. A series of 16‐day forecasts ran every Wednesday at 12Z. For the ensemble configuration, the 60‐day reforecasts were performed every Wednesday at 12Z, starting from 1st February 2017 and through 24 January 2018.
We compared the atmospheric results of the deterministic configuration (Navy‐ESPCDET) to the stand‐alone operational NAVGEM version 1.4. The spectral resolutions of NAVGEM in Navy‐ESPCDET and NAVGEMv1.4 are different (T681 vs. T425, respectively), and, as noted in section 2.1, there are differences in the physics parameterizations between these models. Because of these differences, insight into how the addition of coupling affects the results is difficult to deduce independently, and these comparisons represent differences between Navy‐ESPCDET and the current operational version of NAVGEM. A large amount of analysis was performed in comparing these systems and a very brief summary of some of these results is shown in the following section.
An example of some lower atmospheric verification is shown in Figure 5 examining the 10‐m wind biases for Navy‐ESPCDET and NAVGEMv1.4. For the 10‐m wind speed, the Navy‐ESPCDET has a lower bias in the tropical region, while the biases in the southern extratropics (15°S to 75°S) are higher than NAVGEMv1.4. The differences in the tropical wind biases may be due in part to the coupling as the SST gradients are more realistic compared to the stand‐alone NAVGEM, and wind stress and SSTs are tightly coupled in this region (Chelton, 2001). The spatial maps show a slight improvement in the Navy‐ESPCDET over the stand‐alone system in the western boundary current locations, such as the Gulf Stream and Kuroshio Extension, and these slight improvements are consistent with more realistic SST gradients in Navy‐ESPCDET compared to NAVGEMv1.4. Southern hemispheric winds near the sea ice edge have a greater bias in NAVGEM‐ESPCDET compared to NAVGEMv1.4.
Enlarge this image.5 Figure. (left) Ten‐meter vector wind‐speed biases for (black) NAVGEMv1.4 from FNMOC operations and (red) Navy‐ESPCDET from 1 January 2017 to 31 December 2017. Averages over (top left) 15°S to 15°N, (middle left) 15°S to 75°S, and (bottom left) 15°N to 75°N are shown. The vector wind‐bias is on the y axis and the forecast hour is on the x axis. (right) Spatial maps of biases for (top‐right) FNMOC operational NAVGEMv1.4 and (bottom‐right) Navy‐ESPCDET averaged over the first seven forecast days. In all panels, biases are calculated using ECMWF analyses as verification over the same time frame.
The verification of tropical storms was performed, but lacked robustness due to the limited number of storms in our validation data set. While the sample size does not support an in‐depth analysis, no significant improvement or degradation for the Navy‐ESPCDET against the stand alone NAVGEM v1.4 was found for tropical storm tracks (not shown).
The standard FNMOC global atmosphere scorecard (Hogan, 2014) is an example of a summary which seeks to provide a more comprehensive atmospheric comparison between the Navy‐ESPCDET and FNMOC operational NAVGEMv1.4. The scorecard was constructed by the U.S. Navy Administrative Model Oversight Panel (AMOP), used to evaluate model upgrades for targeted decision makers, and not necessarily model developers concerned about physics, though the information is useful in determining general improvements in a model. The scorecard is defined by a set of specific variables (e.g., wind, air temperature, geopotential height), levels (e.g., 500 hPa, surface, 100 hPa), error metrics (e.g., anomaly correlation, ACOR; root‐mean‐square error, RMSE), regions (e.g., Northern Hemisphere, tropics), verifying references (e.g., radiosonde observations, self‐analysis, buoy observations), and lead times (e.g., 5‐day forecast). For each scorecard entry, we compute pairwise differences between the error metrics of two sets of forecasts over the comparison period, then use a two‐sided t test with a null hypothesis of no difference (i.e., mean pairwise error difference is zero). To account for serial error correlation in subsequent forecasts, we use a lag‐1 autocorrelation correction to compute an effective sample (Wilks, 2011). If the null hypothesis is rejected (and, for some scorecard elements, if the mean difference exceeds a threshold value), then the specified number of points are assigned: positive points for improvement, negative points for degradation. We use a two‐sided t test to evaluate if the forecast errors are distinguishable, rather than sequential one‐sided tests for improvement or degradation. This strategy is consistent with historical practice and decreases the probability that we will erroneously reject the null hypothesis. A neutral change is common for many experiments such as adding an observing system, or minor changes to the model or data assimilation system. Scores ranging from 0–5 are small, 5–10 moderate, and it is rare to see a score larger than 10. The maximum achievable score is 25.
When comparing forecast lead‐times of 4 and 5 days, NAVGEMv1.4 outperforms Navy‐ESPCDET with an overall scorecard difference of 5. However, these degradations are generally small and become even smaller at longer lead times. When considering the same scorecard and extending the lead‐time to 6 days, Navy‐ESPCDET outperforms NAVGEMv1.4 with positive scores in the tropics. These results reflect the design of the Navy‐ESPC physics suite to perform well at longer lead times. The overall findings are that the Navy‐ESPCDET performances are broadly neutral compared with the operational FNMOC NAVGEMv1.4 performance, while providing some significant improvements in the tropics.
Navy‐ESPCDET ocean output is compared to GOFS 3.1 and GOFS3.5. HYCOM in GOFS 3.5 is identical to the version in Navy‐ESPCDET as they both have 1/25° horizontal resolution, 41 vertical hybrid layers and internal tides, but GOFS 3.5 is not two‐way coupled to the atmosphere, i.e., it is a stand‐alone ocean/sea ice prediction system which is forced by operational stand‐alone NAVGEM forecast fields. HYCOM in stand‐alone GOFS 3.1 has 1/12° horizontal resolution and 41 vertical hybrid layers, which is the same as the version in Navy‐ESPCENS. Another difference between Navy‐ESPC and the GOFS systems is that the GOFS systems include flux corrections to account for biases in the NAVGEM fields (Metzger et al., 2013) and biases that develop in SSTs, while Navy‐ESPC does not include flux corrections.
The temperature and salinity structure in the upper 500 m of the globe (spanning 50°S–50°N) is evaluated against independent (unassimilated) profile observations of various types, e.g., Argo profiles (Roemmich & Gilson, 2009). For a given observation, model outputs are sampled at the nearest model grid point and the bias (mean error, model minus observations) and RMSE are computed as a function of forecast length (Figure 6). Verification is further broken down into various depth ranges. All systems show an overall cool bias and RMSE that grow as a function of forecast length. The highest bias and RMSE are in the depth range of the thermocline (50–150 m). The positive impact of atmospheric coupling can be seen in the shallowest depth range (0–50 m) in that Navy‐ESPCDET shows lower bias and RMSE than the stand‐alone ocean forecast systems.
Enlarge this image.6 Figure. Temperature (°C) bias (mean error, model minus observation) (top row) and RMSE (bottom row) evaluated against unassimilated profile observations as a function of forecast length (days) for the global region spanning 50°S–50°N for calendar year 2017. The analyses are grouped by depth ranges: 0–50, 50–150, 150–500, and 0–500 m (first through fourth columns, respectively). The red curves are Navy‐ESPCDET, cyan curves are Navy‐ESPCDET persistence, blue curves are GOFS 3.1, and black curves are GOFS 3.5.
Accurate knowledge of the underwater acoustical environment is also important for naval operations. The three‐dimensional structure of temperature and salinity and the surface mixed‐layer depth (MLD) determine the sound speed profile, which characterizes the acoustical ducts in the ocean. Thus, the forecast system must be able to accurately predict the MLD, sonic layer depth (SLD), and other acoustical proxies. The MLD is defined by a density difference of 0.15 kg/m3 between the surface and a given depth if both temperature and salinity profiles are available or a difference of 0.5°C between the surface and a given depth if only T profiles are available. The SLD is the distance between the surface and the depth of the sound speed maximum, which is often, but not always, at the base of the mixed layer. The MLD and SLD bias and RMSE as a function of forecast length are shown in Figure 7. All systems show a shallow bias with similar performance and all are superior to persistence forecasts of Navy‐ESPCDET analyses, i.e., an unchanging forecast. RMSE grows with increasing forecast length but with similar trends for all systems. However, for this particular metric both stand‐alone GOFS show slightly better performance to Navy‐ESPCDET. The reasons for this are not yet fully understood.
Enlarge this image.7 Figure. Mixed layer depth (m, left column) and sonic layer depth (m, right column) bias (top row) and RMSE (bottom row) against unassimilated profile observations vs. forecast length (days) for Navy‐ESPCDET (red), Navy‐ESPCDET persistence (cyan), GOFS 3.1 (blue), and GOFS 3.5 (black). Note that the results start at 24 hr after the forecast. Hence, differences between the start of the lines.
For sea ice validation, we compared the Navy‐ESPCDET ice output to GOFS 3.1 and GOFS 3.5. A major difference between GOFS3.5 and the other systems is that GOFS 3.5 uses CICEv5.1.2 while GOFS 3.1 and Navy‐ESPC use CICEv4. This newer CICE version includes the correction of known software bugs, two new melt pond parameterizations, an improved parameterization for form drag, use of sea ice velocity in the coupling updates for high frequency coupling, and Elastic‐Anisotrophic‐Plastic rheology to name a few of the changes. The reason for the difference in CICE versions was one of timing, i.e., the ESPC code was locked down for testing before the GOFS 3.5 code.
Figure 8 shows the sea ice edge error compared to the independent analysis of the National Ice Center (NIC) (Partington et al., 2003) as a function of forecast day for the Pan‐Arctic and multiple Northern Hemisphere subregions shown in the inset. GOFS3.5 with CICEv5.1.2 generally outperforms both GOFS 3.1 and Navy‐ESPCDET. However, at forecast lead‐times of 1 and 2 days, Navy‐ESPCDET has similar skill to GOFS 3.5 for the GIN Sea, Barents/Kara Seas, and Beaufort Sea subregions. Note for the Pan‐Arctic and some subregions the error is higher at the nowcast time and decreases as forecast length increases (out to about day 2) before eventually increasing again. We hypothesize this is due to the order in which the Interactive Multisensor Snow and Ice Mapping System (IMS) (Helfrich et al., 2007) sea ice mask is applied, i.e., after the NCODA sea ice analysis completes. IMS typically does not extend as far equatorward as the NIC sea ice edge and so the nowcast time may have higher error. Subsequent testing has shown that bringing the IMS mask inside of NCODA and only applying it to passive microwave sensors that are known to suffer from nondifferentiation of sea ice and summertime melt ponds improves ice edge error performance at the nowcast time.
Enlarge this image.8 Figure. Northern Hemisphere sea ice edge error (km) vs. the independent NIC ice edge as a function of forecast length for (red) Navy‐ESPCDET, (blue) GOFS 3.1, and (black) GOFS 3.5. The geographic locations of the regional analysis are shown in the map at the bottom. Note that Navy‐ESPCDET and GOFS 3.1 use CICEv4.0 whereas GOFS 3.5 uses CICE v5.1.2.
Sea ice drift from the three systems is compared against ice‐bound drifting buoys from the International Arctic Buoy Program (
Enlarge this image.9 Figure. Monthly (top) speed mean absolute error (cm/s), (middle) RMSE (cm/s), and (bottom) direction vector correlation of unassimilated IABP ice‐bound drifting buoys for the Pan Arctic region compared against GOFS 3.1 (blue), GOFS 3.5 (black), and Navy‐ESPCDET (red) for tau 012 at all ice concentrations. The vector correlation is based on Crosby et al. (1993) with 0 indicating no correlation and 2 indicating perfect correlation. This analysis spans all days of 2017.
As noted above, one goal of the NAVGEM physics updates was to obtain skillful long forecasts of the MJO; therefore, we compare Navy‐ESPCENS MJO predictions to state‐of‐the art systems from other NWP centers. In order to evaluate the skill of the MJO forecasts, we use the Real‐time Multivariate MJO (RMM) index (Wheeler & Hendon, 2004). This commonly used index has two components. RMM1 is positive when the MJO is active over the Maritime Continent and negative when the MJO is active over South America to Africa. RMM2 is positive when the MJO is active over the Pacific Ocean and negative when the MJO is active over the Indian Ocean.
We compare Navy‐ESPCENS to the FNMOC operational atmosphere only NAVGEM ensemble transform (NAVGEM‐ET) (McLay et al., 2010; Reynolds et al., 2011), Australian Bureau of Meteorology (BOM), Meteo‐France (CNRM), the European Centre for Medium‐range Weather Forecasts (EC), and the NOAA Climate Forecast System version 2 (CFSv2) (Table 4). We downloaded the forecasts from the other centers from the Subseasonal to Seasonal (S2S) database (Vitart et al., 2017). The models from the other centers vary in ensemble size as noted in Table 4. In addition, we show results for the Navy‐ESPCDET configuration. Verification was derived from NOAA satellite‐derived Outgoing Longwave Radiation, and the ECWMF ERA‐Interim reanalyses.
TableDescription of Modeling Systems Used in the MJO and Teleconnection Intercomparisons| Model acronym | Modeling Center | # of ensemble members | Frequency of forecasts |
| BOM | Australian Bureau of Meteorology | 33 | 3.5 days from 2017‐01‐04 |
| CFSv2 | NOAA Climate Forecast System Version 2 | 16 | Daily from 2016‐12‐31, but match Navy‐ESPCENS |
| CNRM | Meteo‐France | 51 | Weekly from 2017‐02‐01 |
| EC16 | European Centre for Medium Range Weather Forecasting (ECMWF) | 16 | 3.5 days from 2017‐01‐02 |
| EC51 | European Centre for Medium Range Weather Forecasting (ECMWF) | 51 | 3.5 days from 2017‐01‐02 |
| NAVGEM‐D | NAVGEM Operational Deterministic System | 1 | Weekly from 2017‐02‐01 |
| NAVGEM‐ET | Navy NAVGEM Operational Ensemble System | 20 | Weekly from 2017‐02‐01 |
| Navy‐ESPCDET | Deterministic Navy‐ESPC | 1 | Weekly from 2017‐01‐04 |
| Navy‐ESPCENS | Ensemble Navy‐ESPC | 16 | Weekly from 2017‐02‐01 |
Note. Dates are in the format year‐month‐day.
Figure 10 shows violin diagrams indicating the distributions of the day when the ACOR of the individual ensemble member falls below 0.6 for RMM1 and RMM2. The Navy‐ESPCENS ensemble individual members perform exceptionally well in terms of RMM1 ACOR compared to the other centers. In terms of the ensemble mean (''X''s in Figure 10), Navy‐ESPCENS performance is slightly below ECMWF, comparable to CNRM, and better than the other centers and NAVGEM ensemble. Navy‐ESPCDET also performs relatively well.
Enlarge this image.10 Figure. Violin plots showing the distribution of the day that the forecast of RMM1 (top panel) and RMM2 (bottom panel) anomaly correlation (ACOR) falls below 0.6. Each ensemble member is treated as an individual (deterministic) forecast, and the width of the violin is proportional to the number of ensemble members that fall below 0.6 on that day. The probability density is represented by the width of the plot and is proportional to the number of members that have the value. The box‐and‐whisker plot inside the violin plot shows the interquartile range with the whiskers out to the 95% confidence intervals. The ensemble mean is represented by the “X”. Note for RMM2 that the Navy‐ESPCDET forecast is only 16 days and the ACOR did not fall below 0.6 before the end of the forecast.
For RMM2, the performance of individual ensemble members of Navy‐ESPCENS is not quite as good as ECMWF, but on average better than the other centers. For the ensemble mean (“X''s), Navy‐ESPCENS outperforms CFSv2 and BOM. For several of the forecast systems, the average RMM2 ACOR does not fall below 0.6 before the end of their forecast range, hence the lack of “X''s for these systems in the bottom panel of Figure 10.
The MJO forecasts have ACOR skill above 0.6 for about 20 days for both indices, which is comparable to other centers. It should be noted that the sample size here is quite small given that the MJO is not always active and the time period considered is only a year. Thus, these results should be used to give a general idea of performance only. When examining 17 years of Navy ESCP forecasts produced for the SubX experiment, Janiga et al. (2018) and Kim et al. (2019) found that the Navy‐ESPCENS system overestimates the amplitude of the MJO whereas most other models underestimate the amplitude of the MJO. The ACOR of the RMM1 and RMM2 is more indicative of predictions of the phase of the MJO where other metrics such as the RMSE of RMM1 and RMM2 also reflect errors in MJO amplitude.
In addition to the MJO, we anticipated that teleconnection patterns might be predictable on longer time scales compared to fields at one location because of their large scale. Teleconnection patterns evolve on relatively slow time scales and as a result are more predictable than the atmospheric state at any single point. The teleconnection patterns considered include the Arctic Oscillation (AO) and Antarctic Oscillation (AAO), which provide information on the strength and waviness of the northern and southern jet streams, and the Pacific North American Oscillation (PNA) and North Atlantic Oscillation (NAO), which provide information on dominant modes of variability over these respective regions. The models that we compare these indices to are the same as the MJO and listed in Table 4. For MJO, however, the verifying scores were computed using NOAA interpolated OLR (Liebmann & Smith, 1996) and ERA‐Interim (Dee, 2011) zonal and meridional wind fields, while the other teleconnection patterns' verifying scores were computed using the National Centers for Environmental Prediction (NCEP) reanalysis height fields.
Figure 11 shows Violin plots describing the distribution of the forecast day when the ACOR skill score of these ensemble forecasts (considered individually, i.e., in a deterministic sense) drops below 0.6. The individual Navy‐ESPCENS forecasts are comparable or better than the performance of the other systems. The day the Navy‐ESPCENS ensemble mean forecast skill falls below 0.6 is also included for these systems, and is indicated by the ''X''. Using this measure, the Navy‐ESPCENS ensemble mean remains skillful usually about a half‐day or day longer than the median of the individual forecasts, while other centers are usually able to gain 2, sometimes 3, days in skill using the ensemble mean. NAVGEM‐ET gain in skill between deterministic and ensemble for the AO and PNA is particularly large, attesting to the merits in the ensemble design. Thus, the ensemble mean skill between NAVGEM‐ET and Navy‐ESPCENS is comparable as NAVGEM‐ET has a slight edge for the AO, and Navy‐ESPCENS has an edge for the AAO. The skill of the Navy‐ESPCDET is comparable to Navy‐ESPCENS, and generally within a day difference of the ensemble mean score.
Enlarge this image.11 Figure. Violin plots showing the distribution of the day that the forecast of the large‐scale oscillation (as denoted in panels) ACOR falls below 0.6. Each ensemble member is treated as an individual (deterministic) forecast, and the width of the violin is proportional to the number of ensemble members that fall below 0.6 on that day. The probability density is represented by the width of the plot and is proportional to the number of members that have the value. The box‐and‐whisker plot inside the violin plot shows the interquartile range with the whiskers out to the 95% confidence intervals. The day on which the ensemble mean forecast ACOR falls below 0.6 is denoted by the ''X''.
Using the rule‐of‐thumb threshold of 0.6 for useful ACORs, the predictive skill for these patterns typically dips below that threshold between 8 and 10 days. We also wish to note a caveat when interpreting the average skill of the individual members. The relatively good performance of the individual members of Navy‐ESPCENS is in part due to a relative lack of spread in the ensemble. Figure 12 displays the ratio of average variance to average squared error of the ensemble mean, referred to as the spread‐skill ratio, for the Navy‐ESPCENS and NAVGEM ET for selected variables and regions. Here the errors are calculated using ECWMF ERA Interim reanalyses. Neglecting analysis error variance, an ideal ensemble would exhibit a spread‐skill ratio of one. As shown, the Navy‐ESPCENS is less dispersive compared to the NAVGEM ET at lead times before ~15 days in the forecast. This lack of ensemble spread means that the best members in the Navy‐ESPCENS tend not to be as skillful as the best members at other centers, and likewise the worst members in the Navy‐ESPCENS ensemble tend to have smaller errors than the worst members of the other centers. As noted, there is currently a dedicated effort to improve the spread of the Navy‐ESPC ensemble.
Enlarge this image.12 Figure. Average ratio of ensemble variance and squared‐error of the ensemble mean (Spread Skill Ratio) for the (black) NAVGEM ET and the (red) Navy‐ESPCENS. Error is calculated from ERA‐interim analysis. The x axis is the forecast length of the forecast. Four variables for specific regions are shown: (a) 10‐m winds in the North Pole region (65–90°N, 0–360°E), (b) 2‐m air temperatures in the North Pacific Region (30–60°N, 130–240°E), (c) 2‐m air temperature in the Tropical Atlantic Region (20°S–20°N, 290–360°E), and (d) 2‐m air temperature in the Tropical West Pacific Region (20°S–20°N, 110–200°E). The dashed horizontal line at 1 represents a perfect spread‐skill ratio. Note the forecast length is only to 504 hr because this is the forecast length for the NAVGEM ET ensemble.
When diagnosing the 60‐day forecasts of the ocean in Navy‐ESPCENS, the primary reference comparison is climatology. Here we use the Generalized Digital Environmental Model (GDEM) version 4 (Carnes et al., 2010) that is the U.S. Navy's gridded monthly full‐depth climatology of temperature and salinity and their standard deviations on a 0.25° global grid, and was built from 8.3 million profiles spanning multiple decades.
Ensemble spread of Sea Surface Height is first examined to understand the spatial variability created by the ensemble method (Figure 13). Larger spread concentrated in the vicinity of the major current systems associated with the Gulf Stream, Kuroshio (off Japan), and, in the southern hemisphere, the Antarctic Circumpolar Current (ACC), Agulhas Current (off southeast Africa), Eastern Australian Current. The larger spread indicates growing dynamical instabilities associated with the mesoscale features. The ensemble spread grows with forecast length as evidenced from the spread at 240‐hr forecast.
Enlarge this image.13 Figure. The standard deviation of sea surface height (SSH) in meters for the 16 members of Navy‐ESPCENS at forecast hours of (top) 24 and (bottom) 240 hr. Forecasts are described in section 4.
Ocean current validation from Navy‐ESPCENS was performed using independent drifting buoy observations with a drogue at 15 m depth from the NOAA Global Drifter Program (
Enlarge this image.14 Figure. Monthly upper ocean (15 m) speed root‐mean‐square error (RMSE, cm/s) (upper left) and vector correlation (upper right) at the nowcast time and speed RMSE (lower left) and vector correlation (lower right) as a function of forecast length for GOFS 3.1 (red), the control member of Navy‐ESPCENS (blue), and the Navy‐ESPCENS ensemble mean (black) against independent drifting buoy observations. The statistics span 2017.
We show the extended‐range ensemble‐mean temperature bias and RMSE of the top 500 m of the water column in Figure 15. Navy‐ESPCENS mean forecast has skill relative to the climatology, which is indicated by a lower RMSE than the GDEM4 climatology, out to about 35 days in the forecast, as indicated by the crossing of the RMSE lines in the left panel. The RMSE of the control member (without perturbed observations) crosses climatology at a lead‐time of 10 days, which shows the benefit of an ensemble for this metric. The right panel shows the RMSE, bias, and absolute error for the Navy‐ESPCENS ensemble mean and the GDEM4 climatology for the upper ocean to 500 m. Navy‐ESPCENS ensemble mean has a small negative bias at the initial time that decreases overall and becomes positive near the surface, which is consistent with the longwave radiation bias (not shown), although the total heat flux drives the surface temperature. The GDEM4 climatology has a consistent and somewhat larger negative bias. The RMSE of the Navy‐ESPCENS mean at the start of the forecast is noticeably lower than the GDEM4 reference, grows quickly in the upper thermocline, and at 100 m becomes larger than climatology (>1.3°C) between forecast days 26 and 40. Vertical profiles of temperature RMSE and bias averaged over the 60‐day forecast show the overall better performance of the ensemble mean compared to GDEM4. The time average spread to RMSE ratio for temperature in the upper 500 m is 0.60 and increases to 0.67 after subtracting bias from RMSE (RMSE − |ME|), suggesting the Navy‐ESPCENS is accounting for 67% of uncertainty.
Enlarge this image.15 Figure. (left) Mean RMSE of ocean temperature over 8–500 m depth as a function of forecast lead time relative to unassimilated profiles for the ensemble (black line), control member (red line), and GDEM4 climatology (green line), and (RIGHT) temperature bias (ME), RMSE, and mean absolute error (MAE) for the upper 500 m for 1‐ to 60‐day forecast time for the ensemble mean (top row) and the GDEM4 monthly climatology (bottom row).
In situ surface temperature observations, mostly from drifting buoys and ships of opportunity, are used to evaluate the ensemble forecasts along with the GDEM4 climatology (Figure 16). We are aware of larger errors in ship SST observations compared to errors in the drifter observations. However, given the majority of the observations used in our validation come from drifters (~76%), the conclusions drawn are not likely impacted. Note that number of surface observations is significantly larger (~30,000 profiles/day) than profiles observations (~2,500/day). The RMSE of the ensemble mean (black line), control member (unperturbed member) (red line), and GDEM4 climatology are shown in the left panel of Figure 16 and the right panel shows the bias and the standard deviation. For comparison, the errors from 15 ensemble members are shown as gray symbols. Navy‐ESPCENS has lower bias and RMSE than climatology at all lead times out to 60 days. Navy‐ESPCENS bias is near zero at lead times up to 40‐days, which is significantly smaller than GDEM4 at 0.4°C. The control member has skill relative to the climatology out to 20 days in the forecast. The ensemble standard deviation (spread) growth is consistent with the RMSE. However, Navy‐ESPCENS is underdispersive as indicated by the standard deviation lower than RMSE, suggesting that dispersion of the ensembles does not account for the uncertainty defined by climatology. The time‐average spread to RMSE ratio of ~0.72 (ranges between 0.52 and 0.76) suggests that the ensemble is accounting for 72% of uncertainty.
Enlarge this image.16 Figure. (left) RMSE of surface temperature relative to ship and drifting buoy observations for the ensemble mean (back line), control member (red line), and GDEM4 climatology (green line). Gray symbols are RMSE of the 15 ensemble members treated independently. (right) The bias (ME) and standard deviation (STD) as a function of forecast lead time. Dashed lines indicate GDEM4 climatological ME (green) and RMSE (black).
To evaluate the spatial distribution of ensemble temperature forecast skill, the observations and model are binned into boxes spanning 40° of longitude and 20° of latitude for temperature observations between 8 and 500 m and 10° by 10° for SSTs. For each bin, the RMSE is calculated with respect to the ensemble mean (EM), control member (M = 0), and climatology. Next, the RMSE in each bin is averaged over 8–500 m depth or SST. Last, the forecast day on which the ensmble mean RMSE crosses the climatology is identified for each bin, and plotted in Figure 17. We excluded bins with fewer than 50 matchup comparisons for the temperatures from 8 to 500 m (Figure 17, left) and fewer than 100 comparisons for the SSTs with surface data (Figure 17, right). We smooth the time‐series in each bin using a 5‐day boxcar method before finding the day of crossing, and in the event of multiple crossings, the earliest occurrence is selected. The ensemble mean performs better than climatology out to 60 days in the Pacific Ocean except in the Kuroshio Extension east of Japan for both near surface (8 to 500 m) temperatures and SSTs. The western boundary current regions have larger ensemble spread compared to other regions (Figure 13), which may be limiting the predictability. The ensemble mean seems to have a slightly lower predictability to 25–35 days relative to climatology in the Atlantic, Northern Indian Ocean, and Indonesian Seas. The control member has lower predictability than the ensemble mean almost everywhere in the globe (not shown).
Enlarge this image.17 Figure. Spatial distribution of the forecast day when the ensemble mean RMSE crosses the climatological RMSE for ocean temperature in the top 500 m (left) and surface temperature (right) for the Navy‐ESPCENS forecasts.
For the temperature observations from 8 to 500 m, the relatively lower predictability in the western Arabian Sea, Bay of Bengal, and Indonesia regions can be attributed to the monsoon winds. Part of the lower predictability in the Indian Ocean regions likely stems from the dominance of the seasonally reversing monsoon. The gap winds along the west coast of Mexico appear to influence the SST predictability, lowering predictability out to about 15 days. Lower SST predictability is also evident in the Agulhas current region (20–25 days), where the generation of mesoscale eddies dominates the SST variability. When assimilating satellite altimeter data, HYCOM accurately represents Agulhas eddy shedding; however, in forecast mode, it has a tendency to produce and shed more eddies into the South Atlantic than are observed. This is a common unrealistic phenomenon in most eddy‐resolving ocean models (Thoppil et al., 2011). The control member has lower predictability than the ensemble mean almost everywhere in the ocean (not shown).
Extended‐range forecasts of sea ice extent are evaluated using two metrics. The first is the Integrated Ice Edge Error (IIEE) (Goessling et al., 2016), which is the total of the area over which the model has overforecast the extent of the ice edge plus the area over which the model has underforecast the extent of the ice edge. Here we define the ice edge as 15% ice concentration. We also consider a Brier Skill Score (BSS) (Brier, 1950) for ice concentration below 15% following Wayand et al. (2019). The BSS is the degree of improvement of the Brier Score (BS) of the Navy‐ESPCENS forecast over the BS produced using a persisted analysis or climatology. The BSS of one is a perfect forecast, while a BSS of 0 indicates there is no improvement as compared to persistence or climatology.
For verification, we follow Hebert (2015) and use self‐analysis as truth. In addition, the climatological forecasts are produced using SSMR SSM/I data from the National Snow and Ice Data Center (NSIDC) (NSIDC 0051) from the 2007–2017 time‐period. Including analyses from years prior to 2007 may bias the climatology on the high side given the recent observed trends in sea ice. We calculate these scores for both the Arctic and Antarctic, and consider the scores for the entire 2017 period, and for each individual season.
Figure 18 shows the IIEE and BSS for both the Arctic and Antarctic as a function of forecast day averaged for the entire 2017 test period. For IIEE, Navy‐EPSCENS outperforms climatology out to 32 days in the Arctic and out to 40 days in the Antarctic. It outperforms a persistence forecast in the Antarctic for the entire time. In the Arctic, Navy‐EPSCENS outperforms persistence past 5 days, and is comparable to (slightly worse than) persistence during the first 5 days. In terms of BSS, the Navy‐EPSCENS outperforms both persistence and climatology in both the Arctic and Antarctic for the entire 60‐day forecast.
Enlarge this image.18 Figure. (a, b) Integrated Ice Edge Error (106 km2) for the Arctic and Antarctic, respectively, for Navy‐ESPC ensemble mean (maroon), persisted analyses (pink), and 2007–2017 climatology (gray). Range of skill of individual ensemble members indicated by shaded region. (c, d) Brier Skill Score for ice concentration less than 15% for the Arctic and Antarctic, respectively, for Navy‐ESPC ensemble as compared to persistence (pink curve) and 2007–2017 climatology (gray curve). BSS range for individual ensemble members indicated by pink and gray areas.
However, the skill depends on season. For the Arctic, Navy‐EPSCENS retains skill above climatology for the entire 60‐day forecast period for the BS (Table 5). For the Antarctic, Navy‐EPSCENS becomes less skillful than climatology after 46 days for both September‐October‐November and December‐January‐February for the BS. There is considerable seasonal variation in the IIEE metric, ranging from a low of 11 days during June‐July‐August to a high of 46 days for March‐April‐May in the Arctic. In the Antarctic, the Navy‐ESPC becomes less skillful than climatology at 28 days for June‐July‐August and September‐October‐November, but retains skill above climatology for the entire 60‐day forecast in March‐April‐May.
The IIEE metric can be broken down into overpredictions and underpredictions and the largest signal in the results is the tendency of the Navy‐ESPCENS to overpredict sea‐ice extent in the Antarctic for all months except March‐May‐April, and during the other seasons, the overprediction of sea ice is 3 to 5 times larger than the underprediction (not shown). In the Arctic, the overprediction and underprediction areas are generally more balanced, although June‐July‐August have underpredictions that are about 2.5 times as large as overprediction area, which is contributing to the relative high IIEEs during these months.
TableForecast Day for Navy‐ESPCENS Sea Ice Metrics per Season| Statistics for hemisphere | MAM | JJA | SON | DJF | All |
| Arctic | |||||
| BSS < 0 | >60 | >60 | >60 | >60 | >60 |
| IIEE > Climatology | 46 | 11 | 35 | 37 | 32 |
| Antarctic | |||||
| BSS < 0 | >60 | >60 | 46 | 46 | >60 |
| IIEE > Climatology | >60 | 28 | 28 | 39 | 39 |
Note. Brier Skill Score (BSS) < 0 represents the day when the ratio the Navy‐ESPCENS Brier Score over the climatology Brier Score is less than zero. The Brier Score is calculated with respect to 15% sea ice concentration. Integrated Ice Edge Error (IIEE) > Climatology represents when the IIEE of the Navy‐ESPCENS becomes larger than the IIEE from climatology. The top row represents the season in which the forecast started: MAM = March‐April‐May; JJA = June‐July‐August; SON = September‐October‐November; DJF = December‐January‐February; and all months.
Navy‐ESPC is the first global coupled atmosphere‐ocean‐sea ice prediction system developed at NRL. It is a new tool for predicting operational and strategic environmental information for the Department of the Navy. The component models include the atmosphere—NAVGEM; ocean—HYCOM; and sea ice—CICE. Since a motivation for developing Navy‐ESPC is subseasonal forecasting, we introduced physics updates into NAVGEM to improve the subseasonal simulation of equatorial phenomena, particularly the MJO, by updating NAVGEM's deep convection physics from the SAS scheme to a modified Kain‐Fritsch scheme. In addition, we implemented the COARE surface flux algorithm into NAVGEM to provide more consistent air‐sea fluxes with HYCOM, as HYCOM uses the same COARE algorithm. No major physics updates were included in HYCOM or CICE component models to accommodate the longer forecast time horizon. Data assimilation is weakly coupled with NAVDAS‐AR for the atmosphere and NCODA for the ocean/ice components.
In this paper, we described the deterministic (Navy‐ESPCDET) and ensemble (Navy‐ESPCENS) configurations. For Navy‐ESPCDET, we analyzed a database of 16‐day reforecasts spanning January 2017 to December 2017. For Navy‐ESPCENS, we analyzed a 16‐member 60‐day ensemble of reforecasts spanning February 2017 to January 2018 starting every Wednesday at 12Z.
For the Navy‐ESPCDET configuration, the main findings are as follows:
- For the FNMOC global forecast scorecard, operational NAVGEM performs slightly better than Navy‐ESPCDET for 4‐ to 5‐day forecasts. However, the Navy‐ESPCDET outperforms operational NAVGEM at 6‐day forecasts. This reflects the design priorities for subseasonal forecasting in Navy‐ESPC compared to the short‐term (i.e., 5 days) forecasting for NAVGEM.
- For the ocean, Navy‐ESPCDET performance is broadly comparable to the current operational system at tactical forecast lead times (0–7 days). However, performance of a few key acoustical measures (SLD, below layer sound speed gradient) important to Navy operations show a small degradation.
- Sea ice edge errors in Navy‐ESPCDET are generally larger than the sea ice edge errors in GOFS 3.5. This is likely because GOFS 3.5 uses an updated version of CICE (version 5.1.2) compared to CICE version 4 in Navy‐ESPCDET.
For the Navy‐ESPCENS configuration, the main findings are as follows:
- Deterministic forecasts of the MJO are comparable to other leading centers (and better than NOAA CFSv2), with the ACOR dropping below 0.6 level after 20–30 days of the forecast.
- Deterministic forecasts for the AO, AAO, NAO, and PNA are also comparable to other leading centers and have ensemble mean ACORs greater than 0.6 in the 8‐ to 12‐day period depending on index.
- The gain realized by using an ensemble mean over a deterministic forecast is smaller in Navy‐ESPC than in other systems, illustrating the current substantial atmospheric underdispersion in the Navy‐ESPC ensemble system.
- Global ocean temperatures averaged between 8 and 500 m in depth have RMSE values lower than climatology out to about 30 days in the forecast.
- Ensemble mean forecasts of SST compared to surface in situ observations have lower RMSE values than climatology out to 60 days, while deterministic (i.e., single member) forecasts RMSE values are lower than climatology out to 20 days.
- Pan Arctic and Antarctic sea ice extent diagnostics show skill over climatology from 11 to 60 days depending on metric and season. The BSS showed more skillful results compared to the IIEE metric, and Navy‐ESPCENS was more skillful in the northern hemisphere winter than summer for Arctic sea ice extent.
As noted above, a limitation of the Navy‐ESPCENS system is that it is underdispersive, particularly for atmospheric variables. As expected, the EDA system used is not sufficiently accounting for uncertainty in the initial state (i.e., underdispersive), and currently no method to account for model uncertainty is incorporated. Improving ensemble design is a major area of development for the next iteration of Navy‐ESPC. For the atmosphere, methods successfully employed to account for model uncertainty in the stand‐alone NAVGEM ensemble such as Stochastic Kinetic Energy Backscatter (Reynolds et al., 2011; Shutts, 2005) will be tested in Navy‐ESPC. In addition, we are currently testing Analysis Correction‐based Additive Inflation (ACAI) based on the work of Bowler (2017) (Crawford et al., 2020). Further, to improve the Navy‐ESPC ensemble's sampling of initial condition uncertainty we are evaluating the Relaxation To‐Prior Perturbations (RTPP) method (Whitaker & Hamill, 2012; Zhang et al., 2004) as well as the ET method (McLay et al., 2010). In the ocean, spread increases are currently being tested by introducing stochastic perturbations based on Lermusiaux (2006). We expect we will be able to realize even more benefit from the Navy‐ESPCENS forecasts once issues with underdispersion have been addressed (McLay & Shafer, 2016). We are also currently examining how ocean resolution affects spread in the boundary current regions by running an ensemble system at a ¼°.
In addition to ensemble design, multiple other aspects will be added or developed into Navy‐ESPC. We plan on the addition of the wave model, Watch Watch III (Tolman, 1991; Tolman et al., 2002; Tolman & Chalikov, 1996; WW3DG, 2019), and upgrading CICE to a more recent version. Increased vertical and horizontal resolutions and possible physics upgrades will be tested for NAVGEM and HYCOM. In particular, we are testing updates for simulations of the stratosphere in NAVGEM, as a better simulation of the stratosphere may aid in subseasonal forecasting (Scaife et al., 2005). In addition, ocean tides will be added to the ensemble configuration. For data assimilation, we are exploring methods to share covariances across the air‐ocean interface based on the work of Frolov et al. (2016).
We thank the reviewers who provide detailed comments that aided in this manuscript. We thank multiple personnel that aided in Navy‐ESPC model development and discussions of results. This research was sponsored by OPNAV N2N6E and the Office of Naval Research's Earth System Prediction Capability, PE0603207N. This work would not be possible without computer time granted from Department of Defense's High Performance Computing Modernization Program. Please refer to Table 4 for the time period and the calculation of the indices (e.g., NAO, AO, MJO) for data needed to download. An account with ECMWF is needed to download these data and the availability is at no‐cost for research. If you need an account, please follow the “Log In” button. Access to the Navy DSRC may be obtained through a request to the DoD HPCMP (
The S2S database is hosted at ECMWF as an extension of the TIGGE database and can be accessed online (
© 2021. This work is published under http://creativecommons.org/licenses/by-nc-nd/4.0/ (the “License”). Notwithstanding the ProQuest Terms and Conditions, you may use this content in accordance with the terms of the License.
Details
; Metzger, E Joseph 2
; Reynolds, Carolyn A 1
; Ruston, Benjamin 1 ; Rowley, Clark 2
; Ole Martin Smedstad 3
; Ridout, James A 1
; Wallcraft, Alan 4 ; Frolov, Sergey 5
; Hogan, Patrick 6
; Janiga, Matthew A 1 ; Shriver, Jay F 2
; McLay, Justin 1
; Thoppil, Prasad 2 ; Huang, Andrew 7
; Crawford, William 1
; Whitcomb, Timothy 1 ; Bishop, Craig H 8
; Zamudio, Luis 4 ; Phelps, Michael 3 1 Naval Research Laboratory, Marine Meteorology Division, Monterey, CA, USA
2 Naval Research Laboratory, Ocean Sciences Division, Stennis Space Center, MS, USA
3 Perspecta, Stennis Space Center, MS, USA
4 Center for Ocean‐Atmospheric Prediction Studies, Florida State University, Tallahassee, FL, USA
5 Naval Research Laboratory, Marine Meteorology Division, Monterey, CA, USA; Cooperative Institute for Research in Environmental Sciences (CIRES), University of Colorado Boulder, Boulder, CO, USA
6 National Oceanic and Atmospheric Administration, National Centers for Environmental Information, Stennis Space Center, MS, USA
7 Science Applications International Corporation (SAIC), Monterey, CA, USA
8 School of Earth Sciences and ARC Centre of Excellence for Climate Extremes, University of Melbourne, Parkville, Victoria, Australia