気象集誌. 第2輯
Online ISSN : 2186-9057
Print ISSN : 0026-1165
ISSN-L : 0026-1165
Article
第3世代気象庁/気象研究所結合予測システム(JMA/MRI-CPS3)
平原 翔二久保 勇太郎吉田 拓馬小森 拓也千葉 丈太郎髙倉 寿成金濵 貴史関口 亮平越智 健太杉本 裕之足立 恭将石川 一郎藤井 陽介
著者情報
ジャーナル オープンアクセス HTML

2023 年 101 巻 2 号 p. 149-169

詳細
Abstract

A new operational seasonal forecast system, Japan Meteorological Agency (JMA)/Meteorological Research Institute (MRI) Coupled Prediction System (CPS) version 3 (JMA/MRI–CPS3), has been developed. This system represents a major upgrade from the former system, CPS2. CPS3 comprises atmosphere, land, ocean, and sea ice forecast models and the necessary initialization systems for these models. For historical reforecasts, the atmospheric reanalysis dataset JRA-3Q provides initial conditions for the atmosphere and external forcings for land, ocean, and sea ice analysis. In the operational forecast, to initialize the system in near-real time, JMA's operational atmospheric analysis is used in conjunction with JRA-3Q. Next, the land surface model is initialized using an uncoupled free simulation, forced by the atmospheric analysis. The ocean and sea ice models are initialized with the global ocean data assimilation system MOVE-G3, which incorporates a newly developed four-dimensional variational method for temperature, salinity, and sea surface height and a three-dimensional method for sea ice concentration. Compared with the previous system, the CPS3 forecast model components exhibit ∼ 2–4 times higher resolution: the atmosphere and land models are configured with ∼ 55 km horizontal resolution, with 100 vertical atmosphere layers; and the ocean and sea ice models exhibit a resolution of 0.25° × 0.25°, with 60 vertical ocean layers. The physical processes of the atmosphere are greatly refined in CPS3 relative to CPS2, resulting in improved representation of subseasonal to seasonal scale variability, including the eastward propagation of the Madden–Julian Oscillation, winter blocking highs in the North Atlantic, and coupled atmosphere–ocean variability during El Niño–Southern Oscillation events. Our historical reforecast experiment for 1991–2020 suggests that CPS3 demonstrates greater forecast skills than CPS2. The usability of the model output has been improved in CPS3 by reorganizing the operation schedule to provide daily updates of five-member ensemble forecasts.

1. Introduction

Seasonal forecasts provide an outlook of climate conditions three to six months ahead and are used for planning in agriculture and renewable energy production and for preparation for extreme weather events, when the seasonal cycle differs significantly from normal. Weather and climate affect our lives across national borders. For this reason, a number of collaborative frameworks provide consistent forecasts across countries based on objective information sources. One example is the “Regional Climate Outlook Forums”, led by the World Meteorological Organization (2016), where operational numerical forecasting systems provide an objective and scientific basis for forecasts. Demand for seasonal forecasting services increased dramatically in recent years, and a growing need exists for numerical forecasting systems to become more accurate and easier to use.

The main objective of this paper is to describe a new version of a seasonal forecast system, JMA/MRI–CPS3 [Japan Meteorological Agency (JMA)/Meteorological Research Institute (MRI)–Coupled Prediction System (CPS) version 3; hereafter CPS3], which became operational at JMA in February 2022. Predictability of the climate system on subseasonal to seasonal timescales depends largely on interactions among the earth system components, such as the atmosphere and ocean. JMA introduced its first coupled atmosphere–ocean prediction system in the late 1990s (Yoshikawa et al. 2016). Initially, the system specialized in forecasting the El Niño Southern Oscillation (ENSO) and the forecast coverage was limited to the tropics. Subsequent advances in forecasting techniques and a significant increase in computing resources led to the creation of JMA/MRI–CPS version 1 in 2008 (Takaya et al. 2017), which provided global coverage and seasonal forecasts. CPS3, described in this paper, is the third generation of the system.

The primary motivation for the development of CPS3 was to enhance the seasonal forecasting capabilities of the CPS. The secondary goal was to improve the usability of forecast products and to provide forecast updates on a daily basis. The development was also intended to allow seamless expansion of the system so that it could be used for subseasonal forecasts, which would be a major advancement compared with previous systems. This latter development is the most significant change in CPS3 from the previous system, JMA/MRI–CPS2 (Takaya et al. 2018; hereafter CPS2). Seasonal climate and weather forecasts are used in combination, and information should therefore be consistent across these timescales. The importance of this consistency has been widely recognized in recent years, and numerical weather prediction centers accelerated their efforts toward this end (Saha et al. 2014; MacLachlan et al. 2015; Johnson et al. 2019). For CPS3, the first step was to improve the seasonal forecast model so that it could also be used for shorter, subseasonal timescales of two weeks and longer; this required performance and usability improvements to the existing system.

Section 2 begins with an overview of CPS3, followed by a description of the forecast model, initial conditions, initial perturbations, and operational schedule. In order to provide an intercomparison of forecast skill in later sections, particular emphasis will be placed on the changes from CPS2. Section 3 evaluates the forecast performance based on historical reforecast experiments for a 30-year period (1991–2020). Section 4 provides a summary and discusses future issues.

2. System configuration

2.1 Overview

CPS3 is an ensemble forecasting system that uses a coupled atmosphere–land–ocean–sea ice forecast model. The system comprises two parts: forecast model initialization and prediction calculations. Table 1 compares CPS3 with CPS2.

During initialization, CPS3 generates initial conditions for land, lakes, and the ocean and prepares initial perturbations to be implemented in the coupled model. Here, CPS3 relies on externally produced atmospheric analyses: i.e., JRA-3Q (Japanese Reanalysis for Three Quarters of a Century; Kobayashi et al. 2021) and the JMA Global Analysis (GA; Japan Meteorological Agency 2022). These analyses provide the initial atmospheric conditions and external forcings for the land and ocean analyses.

The forecast model comprises GSM (Section 2.2a) and MRI.COM (Section 2.2b), which have both been configured and improved to be appropriate for seasonal forecasting. The model-coupling library, SCUP (Yoshimura and Yukimoto 2008), is used to update air-sea boundary conditions hourly through the exchange of geophysical parameters between the models at the sea surface.

2.2 Forecast model

a. Atmosphere–land surface model

The atmosphere–land surface model in CPS3 is GSM 2003 (Yonehara et al. 2020), which became an operational weather forecast model at JMA in March 2020. The horizontal resolution is set to TL319 (∼ 55 km), and 100 vertical layers are used, giving around double the resolution of CPS2, which used TL159 (∼ 110 km) and 60 vertical layers. The model top was 0.1 hPa in CPS2 and is raised to 0.01 hPa in CPS3. The model dynamics and representations of physical processes incorporate multiple improvements that have been made since GSM 1011 (Japan Meteorological Agency 2013), on which CPS2 was based (Yonehara et al. 2014, 2017, 2018, 2020). To consider layer structure and snowpack coverage, the representation of land surface processes has been reconstructed in CPS3. Next, to better capture diurnal variability, soil temperature and moisture are multi-layered in CPS3. Surface albedo on sea ice has been refined to account for sea ice thickness and snow depth (Hunke and Lipscomb 2010). The scheme from Iwasaki et al. (1989) that was used to represent orographic gravity wave drag in CPS2 is replaced in CPS3 with that from Lott and Miller (1997), which explicitly accounts for low-level flow blocking and orographic gravity waves. Sub-grid-scale turbulent orographic drag is newly considered in CPS3, following Beljaars et al. (2004) (Kanehama and Yamada 2019). The momentum transport effects of nonorographic gravity waves are represented in CPS3 following Scinocca (2003), which improved the reproducibility of stratospheric quasi-biennial oscillations (Kanehama 2012).

In GSM 2003, cumulus convection is parameterized following Arakawa and Schubert (1974), and the set of dynamic and thermodynamic equations is closed using the predictive equation for cloud-base mass flux from Pan and Randall (1998). In CPS3, dissipation timescales for convective kinetic energy are treated separately for shallow and deep cumulus clouds (20 minutes and 40 minutes, respectively). An empirical dependence of the cumulus entrainment rate on altitude and humidity is introduced in CPS3 (Komori et al. 2020; Bechtold et al. 2008), using a minimum threshold taken from Tokioka et al. (1988). These changes increase the relative contribution of shallow cumulus to the total kinetic energy and make the transition from shallow to deep cumulus more continuous. This improves the dry bias for the mid-troposphere and the optically thin cloud bias in the Inter-Tropical Convergence Zone in the Pacific that were seen in CPS2. Observations showed that the air column becomes gradually humidified from the lower to upper troposphere as the convectively active phase of the Madden–Julian Oscillation (MJO) approaches (Thayer-Calder and Randall 2009), and this is captured by CPS3 forecasts. Kawai et al. (2017) provided an elaborate index that describes the conditions of appearance of marine stratocumulus clouds. When this index is sufficiently large, CPS3 weakens vertical mixing to keep a temperature and humidity inversion near the top of the planetary boundary layer, suppressing dry air entrainment from the free atmosphere to generate clouds. A significant melting bias was found for Antarctic sea ice in CPS2, and a related positive bias for downward shortwave radiation fluxes over the Southern Ocean. CPS3 reduced these by improving the representation of super-cooled liquid clouds in the lower troposphere, as these clouds being common in the Southern Ocean region (Kay et al. 2016; Chiba and Komori 2020). Free convective gusts near the sea surface (Godfrey and Beljaars 1991) and deep convective downdrafts are included in the calculation of ocean latent heat release in CPS3, following Redelsperger et al. (2000). In CPS3, we use Zeng and Beljaars (2005) to solve the heat budget in the warm water layer, while allowing changes in the assumed vertical temperature profile, which improves the reproducibility of the diurnal sea surface temperature (SST) cycle.

CPS2 exhibited a dry bias near the land surface over northern hemisphere continents. To address this, we introduced a fractional land ratio to consider subgrid- scale water surfaces in CPS3, shown for Asia in Fig. 1. All water bodies, including rivers, are treated as isolated lakes that are geographically fixed over time in CPS3, and there is no energy or mass transport between adjacent grid cells. Instead, a simple thermodynamic lake scheme is introduced, which predicts lake ice formation and lake temperature variations through water phase changes and heat transfer between water, ice, and snow. In CPS3, several changes were introduced to the representation of atmospheric radiation processes that were used in CPS2. A set of correction schemes from Hogan and Bozzo (2015) and Hogan and Hirahara (2016) are incorporated into the representation of the surface downward shortwave radiation flux to improve the estimated incident net surface radiation at coarse spatio-temporal resolution in CPS3 (radiative fluxes are computed hourly at a resolution of four grid cells).

Fig. 1.

Subgrid-scale land ratio for grid cells in (a) CPS3 and (b) CPS2. The horizontal resolutions of CPS3 and CPS2 are set to TL319 (∼ 55 km) and TL159 (∼ 110 km), respectively. The land ratio is the land area divided by the area of one atmospheric model grid cell.

Further, a monthly climatology is used for ozone concentration in CPS3, as in CPS2, but the climatology in CPS3 has been updated with the 1981–2010 average from the latest MRI-CCM2 reanalysis (Deushi and Shibata 2011). Observed greenhouse gas emissions are used in CPS3 for calculations up to 2016 and are taken from CMIP6 emission scenario SSP2-RCP4.5 (van Vuuren et al. 2011) for later periods. CPS3 uses a three-dimensional monthly aerosol concentration climatology (Yabu et al. 2017) for both reforecasts and operational forecasts and includes an experimental option to evaluate and include the direct radiative effect from volcanic aerosols provided by the user.

Uncertainties in the model physics are calculated using the stochastic physics scheme in Buizza et al. (1999), where physics tendencies are perturbed with space- and time- dependent random noise during model integration; CPS3 continues to use the scheme as it was implemented in CPS2 (Yonehara and Ujiie 2011; Takaya et al. 2018).

b. Ocean and sea ice model

We use MRI.COM (Tsujino et al. 2017), a community ocean model developed at the Meteorological Research Institute, as the ocean and sea ice model in CPS3. CPS3 uses version 4.6 (v4.6) of MRI.COM, as this was the most recent version available at the time of development. MRI.COM uses the Boussinesq approximation to solve the primitive equations using the finite difference method. The horizontal resolution is refined from 1° × 0.3–0.5° longitude–latitude in CPS2 to 0.25° × 0.25° in CPS3. This resolution is sufficient to resolve the first baroclinic Rossby radius for most regions within 30° of the equator (Hallberg 2013) but is too coarse for full resolution at higher latitudes, and it is therefore referred to as an “eddy-permitting” resolution. In CPS3, we use a generalized orthogonal coordinate system in the Arctic (latitudes north of 64°N) on a tripolar grid with singular points in Siberia [64°N, 80°E], Canada [64°N, 100°W], and at the South Pole. CPS3 uses z* vertical levels (Adcroft and Campin 2004), which can accurately capture flow along steep seafloor topography. The number of vertical ocean layers is increased slightly, from 52 in CPS2 to 60 in CPS3, with enhancement primarily in layers deeper than 1000 m. The sea ice model deals with sea ice advection, formation, growth, and melting using five ice-thickness categories. The processes and numerical treatments for the sea ice scheme in CPS3 remain mostly unchanged from CPS2. Further details are available in Tsujino et al. (2017).

Figure 2 compares the SST forecasts for the eastern tropical Pacific for December 22–26, 1999, from CPS2 and CPS3 with an independent satellite-based SST analysis (Merchant et al. 2014). The La Niña conditions that year meant that low SSTs dominated in the equatorial region, and cold water meandered from north to south with Tropical Instability Waves (TIWs). The cold tongue (low SST region) that extends westward from the Galápagos Islands (∼ 90°W on the equator) is enhanced by the coastal upwelling of the eastward equatorial undercurrent close to these islands (Karnauskas 2007). CPS3, with a higher ocean resolution, is able to reproduce these fine-scale SST features more realistically than CPS2. TIWs have been reported to enhance meridional heat exchange across the equator, providing negative feedback to equatorial SST anomalies during ENSO events (Vialard et al. 2001; An 2008; Imada and Kimoto 2012; Graham 2014). In this case, the northward (or southward) flow carries equatorial cold water away from the equator, weakening the amplitude of La Niña. The ability of CPS3 to capture these dynamic effects is likely to mean that the over-development bias for ENSO, which was a critical issue in CPS2, is improved in CPS3.

Fig. 2.

Five-day mean sea surface temperatures (SSTs; °C) for December 22–26, 1999. (a) CCI SST (Merchant et al. 2014) is shown as an independent SST analysis for comparison with the MOVE-G3 and MOVE-G2 analyses used in (b) CPS3 and (c) CPS2, respectively. CPS data are the average for days 11–15 of the forecast for the control member.

2.3 Initial Conditions

a. Initial conditions for the atmosphere and land surface model

JRA-3Q1 provides the initial atmospheric conditions for CPS3 when run in reforecast mode, and GA provides these for operational CPS3 forecasts. There are several differences between JRA-3Q and GA, including system version, resolution, and data cut-off time. JRA-3Q is based on a low-resolution (∼ 40 km) version of the operational global data assimilation system as of December 2018 (Japan Meteorological Agency 2019). The analysis period is extended forward while keeping the version fixed. The GA incorporates developments conducted since then and, as of March 2022, demonstrates a higher resolution of TL959 (∼ 20 km). It will continue to be updated on a regular basis. Although JRA-3Q lags behind real time by about two days, the GA's preliminary analysis provides initial conditions to CPS3 with a delay of only a few hours in exchange for a short cut-off time. Next, we confirmed that these inconsistencies in the initial atmospheric conditions do not critically affect seasonal forecast performance. However, these must be addressed for land surface and ocean analyses because differences in atmospheric forcing accumulate over time and move the mean states away from those in the forecast. For reforecasts, the initial conditions for the land surface are calculated using a free simulation of the stand-alone land surface model of CPS3 itself, using JRA-3Q surface forcing. Only snow cover is used from JRA-3Q. GA is used to calculate the initial conditions for operational forecasts by branching the long-term analysis cycle for one day only. This approach allows us to bring forward the completion time of the simulation while maintaining its historical consistency. Using its own surface simulation avoids forecast “initial shocks” due to discrepancies among the vegetation in the land models implemented in JRA-3Q, GA, and CPS3. Physical parameters in the lake scheme, which are present only in CPS3, are also initialized.

b. Initial conditions for the ocean and sea ice model

In CPS3, initial ocean and sea ice conditions are taken from the global ocean data assimilation system, Multivariate Ocean Variational Estimation/Meteorological Research Institute Community Ocean Model–Global version 3 (MOVE/MRI.COM-G3; hereafter MOVE-G3).

Table 2 shows the major differences between this and the earlier global ocean data assimilation system MOVE/MRI.COM-G2 (MOVE-G2; Toyoda et al. 2013), which was used in CPS2. In MOVE-G2 and MOVE-G3, gridded SST analyses are assimilated as though they were observation data. MOVE-G3 uses MGD SST (Kurihara et al. 2006), a quarter-degree resolution analysis that includes data from satellite observations, whereas MOVE-G2 uses COBE-SST (Ishii et al. 2005), a one-degree resolution analysis that is based on in-situ observations. In addition to the change to the SST products, the assimilation scheme changed significantly between MOVE-G2 and MOVE-G3. A 4D-Var method is used in MOVE-G3, which deals with inhomogeneous observation times better than the three-dimensional method (3D-Var) used in MOVE-G2 and generates dynamically balanced initial conditions for the forecast model. One issue with using 4D-Var in CPS3 is the high computational cost. This, combined with the higher resolution of the ocean model that must be initialized for CPS3 (relative to CPS2), means that much higher computational resources are required for CPS3. To address this, we perform 4D-Var on a lower resolution grid of 1° × 0.3–0.5° (G3A) and downscale the analyzed fields onto a 0.25° × 0.25° grid (G3F) using Incremental Analysis Updates (IAU; Bloom et al. 1996). This twostep approach is based on Usui et al. (2015) and improves the accuracy of analysis without requiring the computational resources needed for a full-resolution 4D-Var. The design of G3A and G3F is similar to the inner- and outer- models used in incremental 4D-Var schemes, but forecast fields are not passed from the high-resolution model (G3F) to the low-resolution model (G3A) for the first guess, and temperature and salinity fields in G3F are nudged to the G3A analysis fields using the IAU instead of applying the analysis increments of G3A to G3F directly. This results in the additional benefit of avoiding computational instabilities that could arise owing to unnatural currents attributable to inconsistencies in ocean topography between the different resolutions.

Figure 3 assesses and compares the water temperature field in the new and old ocean analysis. Root Mean Square Error (RMSE) is estimated through comparison with ARGO float data that are withheld from the analysis. The comparison suggests that G3A provides better estimates in many regions. In particular, temperatures at 1 m depth show clear improvement (Fig. 3b), which can be attributed to the use of MGD SST and possibly to the introduction of 4D-Var. Close to the sea ice, a checkerboard pattern of large RMSE differences is found, which may be due to the small number of ARGO floats available in that region. The large RMSE in some coastal areas is likely due to the inability of low-resolution 4D-Var to represent coastal upwelling, coastal currents, and the associated nonlinearities. For 100 m depth, improvements are modest in most areas (Fig. 3d). The 4D-Var analysis exhibits considerably more variability than the analysis that uses 3D-Var, which may account for the improvements in RMSE being smaller than expected. Some improvements are found around subtropical gyres in the South Indian and South Atlantic Oceans. As the water temperature at 1 m depth (Fig. 3b) is consistently more accurate in these areas when 4D-Var is used, the combination of improved atmospheric forcing and a better analysis method may have improved the representation of the large-scale circulation.

Fig. 3.

Root mean square error (RMSE) for water temperatures in MOVE-G3 and the percentage difference in RMSE relative to the older system, MOVE-G2, at depths of (a, b) 1 m and (c, d) 100 m. A reanalysis experiment was conducted for the period 2005–2014 using the old and new systems, with data from 20 % of the Argo floats withheld from assimilation into the reanalysis and used to evaluate the RMSE of the water temperature in the reanalysis.

Another major improvement is the introduction of sea ice assimilation (Toyoda et al. 2016). MOVE-G2 did not assimilate any sea ice observations. Instead, it was constrained through the dynamics and thermodynamics of the forecast model that assimilated other observations. One example that exposes the shortcoming of this approach is that the sea ice modeled around Antarctica was underrepresented in response to the positive bias in the incoming shortwave radiation flux in JRA-55 (Kobayashi et al. 2015), and to allow the atmospheric forcing to match the satellite observations, the analysis scheme was required to adopt an empirical bias correction. Also, with the introduction of the improved atmospheric forcing of JRA-3Q and data assimilation of sea ice concentration, MOVE-G3 no longer needs to apply such an empirical correction method. In MOVE-G3, a daily, quarter-degree operational sea ice concentration analysis (Matsumoto et al. 2006) is assimilated using 3D-Var. Although sea ice concentration is the only constraint used in this analysis, other parameters, such as ice thickness and sea surface salinity and temperature in ice-covered regions, are updated in the analysis cycle through forward-model integration. Here, the 3D-Var and IAU for sea ice are performed independently for G3A and G3F so that the coastal topography is represented optimally in both.

To illustrate the impacts of the newly introduced sea ice assimilation, Fig. 4 compares the CPS reforecasts with the climatological sea ice extent for the first forecast month in the Arctic Ocean (See Section 3.1 for details on this reforecast). The agreement of the mean ice edge locations in the analysis and in the CPS products indicates that no unnatural initial drift occurs in the CPS and that the observations are properly assimilated. Figure 4 shows data for September and March because these are the months when Arctic sea ice reaches its smallest and greatest extent, on average. The comparison shows that the bias in CPS2—which does not assimilate sea ice concentration—is very small. This may be because the sub-zero SSTs and sea surface fluxes used for data assimilation effectively controlled the formation and disappearance of sea ice. However, the comparison also shows that assimilating sea ice concentration results in a closer agreement with observations. The impact is particularly clear in September, when the sea ice starts to retreat from the provided initial conditions toward the pole after the CPS2 forecast begins. The positional bias of the sea ice edge is particularly improved in the Greenland Sea (B in the figure) and the Chukchi Sea (C) in CPS3, relative to CPS2. We separately confirm that the improvement in the mean error leads to better anomaly-correlation scores for sea ice concentration itself and for 2 m air temperature. The impacts during the freezing season are relatively small, probably because sea ice formation depends more on the chaotic temporal evolution of the sea surface wind and heat fluxes than on the initial conditions in the model. The improvements are not sufficient to address the under-estimations of sea ice extent in the Labrador Sea (A) and the Sea of Okhotsk (D), although the ice edges agree well with the assimilated observations at the initial state.

Fig. 4.

Climatological ice cover in March and September. Contours show a climatological sea ice concentration of 0.15 in March (blue line) and September (black line) in COBE-SST (analysis). Shading indicates a concentration of ≥ 0.15 in March (light blue) and September (gray) in reforecasts from CPS3 (left) and CPS2 (right). The model climatology for March (September) is constructed from a 10-member ensemble forecast initialized on February 10 and 25 (August 14 and 29, respectively). Regional pattern correlation coefficients for sea ice concentration defined over the Arctic Ocean (40–90°N) are displayed in the lower left of each panel. Capital letters denote the positions of the Labrador Sea (A), Greenland Sea (B), Chukchi Sea (C), and Sea of Okhotsk (D).

To initialize the ocean model once every five days, MOVE-G2 performs a preliminary analysis. In contrast, MOVE-G3 is designed to produce an analysis every day by running five analysis streams and executing preliminary analysis for one of the streams each day, using observations from the last five days. Also, MOVE-G3 implements a “delayed-mode analysis” that waits for observations up to nine days. Both analyses precede each assimilation run by five days. MOVE-G3 uses JRA-3Q and GA data for the atmospheric forcing. The idea behind this is the same as the basis for the aforementioned land model initialization: JRA-3Q is used for delayed analysis, for historical consistency, whereas GA is used for the preliminary analysis because of its immediate availability and consistency with the atmospheric initial conditions. The surface heat flux bias has been reported to be much lower in JRA-3Q than in JRA-55 (Kobayashi et al. 2021), making it a more suitable data source for the atmospheric forcing for the ocean analysis.

2.4 Initial Perturbations

a. Initial perturbations for the atmosphere model

Next, to reflect uncertainties in the atmospheric analysis, small perturbations are added to the initial conditions for the ensemble forecast. CPS3 uses the Breeding of Growing Mode (BGM; Toth and Kalnay 1993; Chikamoto et al. 2007) method to extract a set of fastest growing error modes. For this purpose, the atmosphere-only forecasts are calculated for 24 hours. The norm is defined from the root mean square of the variability of the 500 hPa geopotential height, averaged separately over the northern (20–90°N) and southern (20–90°S) hemispheres and from the 200 hPa velocity potential for the tropics (20°S–20°N). The estimated perturbation patterns are rescaled with positive and negative coefficients and added to the analysis. The rescaling factors are fixed in both the reforecast and operational forecast at 14.5 % of the climatological variability for the 500 hPa geopotential height and at 20 % of the climatological variability for 200 hPa velocity potential; this assumption is made for simplicity. Ideally, the size of the initial spread should change, as the accuracy of the atmospheric analysis is not constant over time. In fact, experiments for the summer 2020 (June–July–August) period show that the spread-skill ratio of the 500 hPa geopotential height for the northern hemisphere is above 2 until the 72nd hour of the forecast, when it should ideally be 1, suggesting that the initial perturbations are too large for the accurate initial conditions of recent years. Although not critical to seasonal forecasting applications, this issue will be addressed in future work.

b. Initial perturbations for the ocean model

To represent uncertainties in the ocean initial conditions, CPS3 uses perturbations that approximate the analysis error covariance structure in the 4D-Var (Fujii et al. 2022). To minimize a cost function, G3A employs a quasi-Newton method, where control variables are iteratively updated (Fujii and Kamachi 2003; Fujii 2005). The size of the updates applied to the control variables and the gradient of the cost function can be used to obtain approximate estimates of the eigenvalues and eigenvectors for the error covariance matrix for the analysis (Niwa and Fujii 2020). In CPS3, initial perturbations are created by combining the estimated eigenvectors after scaling so that their amplitude equals half the analysis increment for the specific day.

Figure 5 shows an example of the spread of the initial ocean perturbations. Compared with CPS2, CPS3 exhibits stronger perturbations that spread over multiple vertical layers, rather than only near the thermocline. The larger spread is due to the arbitrarily chosen scaling factor mentioned above, but the change in the pattern reflects a fundamental improvement in the way that the perturbations are generated. In CPS2, ocean perturbations were created solely from the atmospheric forcing (Takaya et al. 2018). The atmospheric forcing used for this was also used to calculate the atmospheric perturbations through the BGM; therefore, it was not an ideal basis for calculating appropriate ocean perturbations, particularly when the MJO is weak or displaced. Part of the under-representation of the Central Pacific spread in MOVE-G2 in Fig. 5 is due to the fact that the convectively active phase of the MJO was in the Atlantic to Indian Ocean near the end of May 2012 and not in a location that would result in strong perturbations in the Pacific. CPS3 is not affected by these issues and provides a straightforward representation of the uncertainties that are inherent to the ocean analysis.

Fig. 5.

Initial ensemble spread for water temperature at the equator on May 31, 2012. Shading indicates the standard deviation for water temperature perturbations (°C) for (a) CPS3 and (b) CPS2 at the equator on May 31, 2012. The black line indicates the 20°C isotherm, which serves as a guide for the tropical thermocline. The vertical axis represents the depth below the sea surface.

2.5 Operational Schedule

CPS2 performed operational model integrations for up to seven months at a time. This provided the basis for operational ENSO forecasts for the next six months using the Lagged Average Forecast (LAF) method (Hoffman and Kalnay 1983). This configuration was carried over to CPS3, although the operational schedule for CPS3 differs significantly from that of CPS2; in CPS2, 13-member ensemble forecasts were produced every five days. Model integrations started two days after the forecast initial time and completed at three more days later (Fig. 6). In contrast, CPS3 calculates 5-member ensemble forecasts on the same day as the initial forecast time. As described in Section 2.3, the analysis cycle was revised to accomplish this change; using GA data, the initial forecast conditions are calculated with a delay of less than 6 hours (the delay was about 54 hours for JRA-55 and 60 hours for MOVE-G2).

Fig. 6.

Operational schedule for CPS2 and CPS3. Numbers indicate the number of members per initial day, and arrows indicate when to start and complete the model integration. To highlight the differences between CPS2 and CPS3, only the numbers and arrows for forecast initial dates of December 22 and December 27 at 00 UTC are drawn in red and blue, respectively.

This change means that users now receive access to 25 ensemble members for the same five-day period, nearly double the number that were previously available. If no requirement for the five-day interval exists, then the ensemble size and the length of the LAF can be optimized on a daily basis.

3. Verification of CPS3, based on a 1991 - 2020 reforecast

3.1 Reforecast settings

In this section, we briefly compare the forecast skill of CPS3 and CPS2, based on a reforecast for 1991–2020. The same experiment design is used for CPS2 and CPS3 to calculate 5-member ensemble reforecasts that each start at 00 UTC in the middle of the month and at the end of the month (Takaya et al. 2018). The start dates are January 16 and 31, February 10 and 25, March 12 and 27, April 11 and 26, May 16 and 31, June 15 and 30, July 15 and 30, August 14 and 29, September 13 and 28, October 13 and 28, November 12 and 27, and December 12 and 27. These dates are partly determined by the fact that initial ocean conditions were only available for CPS2 once every five days; to facilitate comparison, the same experiment design is applied to the reforecasts from CPS2 and CPS3. The oldest initial date for the reforecasts was set to be 15 days behind the latest initial date for each month, following the operational LAF configuration of CPS2.

For verification of the subseasonal forecasts (Section 3.2), five ensemble members for each initial date are used from Day 1. For verification of the seasonal forecast (Section 3.3), 10 ensemble members from each month are aggregated and used to set the monthly and ensemble means from the beginning of the next month. For example, the forecasts that begin on December 12 and 27 are used to calculate the monthly and ensemble averages of January (Month 0), February (Month 1), and so on. The forecast performances are assessed through comparison with data from JRA- 3Q, MGD SST, NOAA outgoing longwave radiation (OLR; Liebmann and Smith 1996), and GPCP v2.3 (Adler et al. 2018). MGD SST and JRA-3Q data are used for the initialization of CPS3, and this may unfairly benefit CPS3 in the comparison. Therefore, we replaced these data with data from COBE-SST and JRA-55, which are used in CPS2, and found that this made no significant difference to our conclusions from the comparison.

3.2 Subseasonal forecast

The MJO exhibits pronounced subseasonal variability in which active tropical convections travel eastward with an average period of 30–60 days and affect mid- and high-latitude variability through atmospheric teleconnections. Therefore, when evaluating global model performance, beginning with the MJO is natural. Figure 7 shows the longitude-time composite for OLR anomalies, predicted from initial dates when the convectively active phase of the MJO was in the eastern Indian Ocean. The comparison suggests that CPS3 represents the eastward propagation of the active/inactive convection phases well, whereas CPS2 exhibits a bias where the convectively active phase becomes stuck in the western Indian Ocean. Separately, we confirmed that this bias was more strongly observed in boreal summer than in winter. Also, the bias was seen for other initial forecast conditions, so that the convections tended to stagnate once the active phase entered the Indian Ocean. Therefore, the improvement is likely to be due to updated physical parameterizations, rather than to the improved representation of the initial conditions. To confirm the forecast skill, we calculated correlation coefficients for an all-season MJO index (Wheeler and Hendon 2004), which remained above 0.5 until Day 21 for CPS2 and Day 27 for CPS3 (not shown). This score compares favorably with recent numerical forecasting systems (Vitart 2017).

Fig. 7.

Longitude–time composite for outgoing longwave radiation (OLR) anomalies for forecasts starting from a convectively active phase of MJO in the eastern Indian Ocean (Phase-3). The horizontal axis is longitude, and the vertical axis is forecast lead-time [days]. Anomalies are defined as the deviations from the 1991–2020 average. The initial phase is detected according to the definition of Wheeler and Hendon (2004). All initial months were included.

Another improvement in CPS3 relative to CPS2 can be seen in the representation of winter blocking highs in the northern hemisphere (Fig. 8). Both models tend to underestimate the frequency of the blocking highs, and, although this bias is greatly reduced in the North Atlantic in CPS3, little or no improvement is found over the north Pacific in CPS3 relative to CPS2. One reason for the change may be the increased resolution. Previous studies reported that increasing the horizontal and vertical resolutions of atmosphere models results in a better representation of dynamic feedbacks between blocking highs, transient eddies, and the terrain effects of steep mountains, but previous researches showed that such effects are only visible in the Atlantic and that differences are found among models (Anstey et al. 2013; Berckmans et al. 2013; Schiemann et al. 2017). This is consistent with our result. To isolate the impact of the increased resolution in CPS3, we compared our atmosphere model (TL319) with a lower resolution configuration (TL159) and confirmed that the high resolution resulted in a consistent improvement. However, the possibility exists that other changes between CPS2 and CPS3 contributed more to the improvements in the results from the newer model system than the resolution changes, as numerous changes exist that could also affect the representation of blocking highs, such as the gravity wave stress scheme that was introduced in CPS3 (Pithan et al. 2016).

Fig. 8.

Climatological frequency (per day) of blocking highs in (a) JRA-3Q, and differences between JRA-3Q and (b) CPS3 and (c) CPS2. Blocking highs are detected following Scherrer et al. (2006) using the 7-day mean analysis for November–February 1991–2020 and corresponding forecasts, where the central day of the 7-day average belongs to this season at days 4–27. Black lines indicate climatological frequencies with values above 0.05, at 0.05 intervals. Values north of 75°N are not defined in this blocking index.

Next, we compare the anomaly correlation coefficients (ACCs) for weekly averaged 500 hPa geopotential height as a measure of the subseasonal forecast skill (Fig. 9). The scores show a significant improvement in CPS3 relative to CPS2 for weeks 1 and 2 in both hemispheres. Subsequent weeks are omitted, but the significant improvement continues until ACC approaches its lowest value in weeks 3 to 4. For short lead times, the improvements may be partly attributable to the use of the latest reanalysis, JRA-3Q, as initial conditions. For later lead times, the possibility also exists that refinements in the model physics may have contributed further. In the Northern Hemisphere winter season (December–January–February), a peak of score improvement can be found over the North Atlantic (figure not shown). This fact is consistent with the improved climate reproducibility of the blockings in CPS3 (Fig. 8). In addition, the improved MJO (Fig. 7) may also have contributed to the overall mid-latitude scores through remote influences. Kubo and Ochi (2022) compared CPS3 with the operational subseasonal forecast system—the JMA Global Ensemble Prediction System (Yamaguchi et al. 2021)—and reported that CPS3 demonstrates a comparable skill when the same ensemble configuration is used for both.

Fig. 9.

Anomaly correlation coefficients (ACCs) for the weekly mean 500-hPa geopotential height. The ACCs are based on the statistics for forecast lead weeks 1 and 2 in December–January–February or June–July–August 1991–2020. The ACCs between the forecasts and JRA-3Q are calculated for each 2.5° × 2.5° grid and averaged over the northern (20–90°N) or southern (20–90°S) hemispheres. Error bars represent 95 % confidence intervals estimated over 1000 bootstrap trials for all forecast initial dates in each season.

Increasing the ensemble size would make CPS3 more appropriate for subseasonal forecasting; however, this presents a challenge. A forecast ensemble gives the probability of occurrence for future weather and climate conditions with a limited number of samples. Increasing the ensemble size is one way of reducing the estimation error for the probability distribution function. In the LAF approach, the shortfall in the number can be compensated for by including data from forecasts with older initial dates in the ensemble. For slowly time-evolving phenomena such as ENSO, Trenary et al. (2018) reported that greater performance can be expected by extending the LAF length beyond a few days and aggregating more ensemble members. However, the LAF length should be much shorter for shorter lead-time forecasts. Figure 9 can be viewed as a comparison of scores among ensemble members within a 1-week LAF. The ACCs are as high as 0.9 for Week 1 and decrease rapidly to about 0.5–0.7 when the initial date becomes a week older. If we combine the forecasts from these initial dates to form a Week 1 forecast, the forecast skill will accordingly deteriorate considerably. For some applications, even a delay of a few days may be critical. By reducing the delay to the initial forecast time and increasing the effective ensemble size (Fig. 6), CPS3 improved usability on subseasonal timescales compared with CPS2, but further enhancements are needed to make it suitable for wider use in the future.

3.3 Seasonal forecast

ENSO is a major source of predictability on seasonal timescales. Niño indices, defined as the regionally averaged SST in the tropical Pacific, are useful measures that provide a brief overview of a system's seasonal forecasting capability. Figure 10 compares analyzed SST anomalies averaged over NINO3.4 [170–120°W, 5°S–5°N] with those predicted from June for the latter half of the year. The figure shows that CPS2 tends to overdevelop initial ENSO signals, particularly in the mid-2010s, whereas CPS3 tends to avoid this monotonous time evolution. Also, CPS3 represents clearer case-to-case variability in the forecast spread, although this is still not adequate to capture the observations. This indicates that CPS3 can simulate a wider variety of ENSO development scenarios. In Section 2.4, we showed that CPS3 can generate effective ocean initial perturbations, using the June 2012 case as an example. A marked contrast is found between CPS2 and CPS3 for ENSO forecasts that are calculated from this initial month (Fig. 10). Although the newly developed perturbation in CPS3 was not designed specifically to capture ENSO, the possibility exists that it helped to diverge the initial development of ENSO, at least for this case. Another contributing factor is the improved forecast model. As shown in Fig. 7, the bias that caused the MJO to stay at a particular longitude in CPS2 is reduced considerably in CPS3. This allows for a more chaotic time evolution of sea surface wind, which has been reported as a key to successful ENSO predictions (Moore and Kleeman 1999; Kessler and Kleeman 2000). In terms of negative feedback processes, TIWs are represented in more detail in CPS3 than in CPS2 (Fig. 2), thereby suppressing equatorial SST anomalies and the associated diffusion effect of equatorial anomalies. Also, improvements in the shallow cumulus cloud scheme contribute to the suppression of excessive ENSO through negative cloud-radiative feedbacks (Wood and Bretherton 2006), as reported in Chiba and Kawai (2021).

Fig. 10.

NINO3.4 SST anomaly time series for analysis and forecasts. The figure shows forecasts from (a) CPS3 and (b) CPS2 for lead times of 0–6 months (June–December) from June initial conditions (10-member LAFs from May 16 and 31) for each year. The black line is MGD SST, red lines are individual ensemble members, and the blue line is the ensemble average.

To overview the changes in ENSO prediction scores, Fig. 11 compares the ACCs calculated using all forecast cases. It can be seen that CPS3 shows a consistent increase in scores from Month 0 to Month 6, with the difference from CPS2 being statistically significant in the first few months. Previous studies reported that ENSO forecast skill sharply declines in boreal spring (referred to as the “Spring Predictability Barrier”; Webster and Young 1992; Tang et al. 2018). When broken down by forecast initial month, early spring to early summer months (February to June) show a clear improvement in forecasts for summer and later seasons, indicating that the skill decline, or the predictability barrier, appears more slowly in CPS3. The RMSE is significantly lower for all lead times in CPS3 than in CPS2. The change in the forecast spread is small; however, the large reduction in RMSE means that the performance improved from CPS2 to CPS3 in terms of the spread-to-skill ratio. In particular, forecasts initialized in Month 0 demonstrate a significantly smaller RMSE and a larger spread, bringing the spread-skill ratio closer to one. This change is likely to be a result of the improved initial ocean conditions (Fig. 3) and the initial perturbations that were introduced in CPS3. In another comparative experiment, we found that the introduction of 4D-Var to the ocean analysis in CPS3, replacing 3D-Var in CPS2, significantly reduces RMSEs for NINO3.4, especially in forecasts with early lead times up to Month 2. The larger forecast spread for early months may be a result of the perturbation of sub-surface ocean layers in CPS3, which was not implemented in CPS2 (Fig. 5). However, also note that the change in the forecast spread is, on average, smaller at longer lead times. In other words, even though CPS3 gives stronger perturbations near the thermocline, in some cases, the perturbations do not develop as much as expected. Scaling the initial ocean perturbations larger might seem to improve the situation. As it turns out, our model tends to dissipate rather than develop such overly large initial perturbations. Considering that these perturbations are based solely on the errors in the ocean analysis, this may be because the unstable modes induced by the added perturbations do not always match the modes that should develop in the coupled atmosphere–ocean system. The ability to represent errors in initial conditions is in itself a steady progress, but further understanding of the error development process is also needed.

Fig. 11.

(a) Anomaly correlation coefficients for NINO3.4 SST between MGD SST and CPS3 (red) and CPS2 (black). (b) Root mean square error in Kelvin (solid line) and forecast spread (dashed line). Statistics are based on 360 cases extracted from the two initial dates of each month in 1991–2020. Lines show the means of 1000 bootstrap trials and error bars show the 95 % confidence intervals.

Next, to compare the typical pattern of atmosphere and ocean variability during ENSO events, Fig. 12 shows the distribution of SST, precipitation, and sea level pressure regressed onto NINO3.4 SST for December–January–February. The analysis shows that, in this reforecast period, ENSO fluctuations tend to occur mainly in the central Pacific (∼ 150°W). The variability in CPS2 is skewed toward the eastern Pacific, and the forecast SST and precipitation anomalies are stronger than the observations. Also, these biases are still evident in the CPS3 forecasts, but they are improved relative to CPS2. Equatorial SST anomalies are meridionally broader near 150°W in CPS3 than in CPS2, suggesting the greater diffusion effect of SSTs (Vialard et al. 2001; An 2008; Imada and Kimoto 2012; Graham 2014) in the higher-resolution model (Fig. 2). La Niña events are often weaker than El Niño events. This asymmetry is particularly pronounced in the eastern equatorial Pacific. We separately confirmed with the frequency distribution of SST anomalies in NINO3 [150–90°W, 5°S–5°N] that CPS3 represents such nonlinearities more closely to observations. Also, CPS3 reproduces sea level pressure anomalies well in the Indian Ocean (60–120°E) and the western tropical Pacific (120–150°E). Further, we confirmed that lower tropospheric circulation fields, such as 850 hPa winds, are also improved in CPS3 relative to CPS2 over a wide area from the Indian Ocean to the western Pacific, suggesting that the interannual variability of the winter Asian monsoon improved through the improved atmospheric response to ENSO in CPS3.

Fig. 12.

Regression coefficients between NINO3.4 area-averaged and global SST (shading; K K−1), precipitation (black line; mm day−1 K−1), and sea-level pressure (blue line with hatching in the areas above 1.2; hPa K−1) during boreal winter (December–February). The regression coefficients are based on the statistics of the November initial conditions in 1991–2020.

Surface air temperature is one of the primary variables of interest in seasonal forecasts. Figure 13 summarizes the regionally and three-month averaged ACCs for 2 m temperature and related variables. In general, the forecast skill for surface air temperature is unchanged or improved from CPS2 to CPS3 for most regions and seasons in the tropics and the northern and southern hemispheres; however, the error bars are large. Precipitation in the tropics is significantly improved for all seasons in CPS3 relative to CPS2. A separate geographical comparison confirms that the regions where large differences exist in the precipitation ACCs for CPS2 and CPS3 coincide well with areas of high interannual variability around the ITCZ, including the western and eastern tropical Pacific. As in the case of ENSO in Fig. 12, this suggests that areas of atmospheric convection tend to be more accurately located in CPS3 than in CSP2. Precipitation anomalies are associated with circulation responses such as vorticity generation in the lower/upper troposphere and remote influences on the mid-latitudes through atmospheric teleconnections. The 850 hPa stream function consistently shows an overall improvement in the tropics and extratropics from CPS2 to CPS3. In regions where the lower-tropospheric circulation influences surface air temperature variability, such as the northwestern edge of the northern hemisphere subtropical anticyclone, changes in the representation of tropical air–sea circulation patterns may contribute to the improved 2 m temperature scores in CPS3 relative to CPS2. Future work will involve investigating regional differences in scores to better understand the factors that influence them by using case studies during the reforecast period.

Fig. 13.

Anomaly correlation coefficients (ACCs) for the average of months 1–3. The vertical and horizontal axes represent the ACC and variable name and region, respectively. The variable names “Tsurf ”, “Prec”, and “PSI850” denote 2-m air temperature, precipitation, and the 850 hPa stream function, respectively. The region names “TR”, “NH”, and “SH” indicate that ACCs are averaged over the tropics (20°S–20°N) and the northern (20–90°N) and southern (20–90°S) hemispheres, respectively. ACCs are computed for Months 1–3 averages for forecasts starting from each initial month and are summarized by the season to which Month 2 belongs. The error bars represent 95 % confidence intervals estimated over 1000 bootstrap trials for all forecast initial dates in each season.

4. Summary and conclusions

We described CPS3, a new operational seasonal forecast system. The latest atmospheric reanalysis, JRA-3Q, and ocean analysis that incorporates 4D-Var and sea ice data assimilation schemes improved the initial conditions used in CPS3 forecasts relative to those used in CPS2 forecasts. The introduction of a high-resolution forecast model and the refinement of physical processes within the model result in an improved representation of interannual variability over a wide range of timescales relative to CPS2, including for the MJO and for blocking highs and ENSO. These improvements were confirmed by forecast scores in reforecast experiments for 1991–2020. To achieve a large ensemble, the operational configuration for CPS3 continues to follow the LAF method but nearly doubles the effective ensemble size relative to CPS2. In addition, the forecast update interval changed from five days in CPS2 to daily in CPS3. This gives users greater flexibility for configuring their LAF ensemble. This usability improvement is not captured in the forecast scores.

The key objective of this report was to present the basic specifications for CPS3 and to describe the differences between CPS3 and the previous system, CPS2. Only elements considered relevant to the headline scores of the operational subseasonal and seasonal forecasts were included in our evaluation. The seasonal characteristics of atmosphere and ocean circulation in the model were outside the scope of this paper, and future work should provide a detailed analysis of these.

CPS3 was developed primarily for seasonal forecasting applications. To make it more suitable for shorter timescales, using a larger ensemble than the current five members per day would be advantageous. This would allow for a shorter LAF length and would mitigate the deterioration of forecast skill over time. Also, a larger ensemble size is needed to allow reforecasts to more accurately capture past ENSO events (Doi et al. 2019) and to estimate a more accurate model climatology.

We reported that the accuracy of the ocean analysis is compromised in some areas by the insufficient resolution of MOVE-G3. This is mainly due to our adoption of the two-step approach described in Section 2.3b. Although another cost-effective alternative could be explored, a straightforward solution is to simply increase the resolution of the analysis when more computing resources become available in the future. Also, increasing the forecast model resolution to become eddy-resolving (∼ 0.1°) would be beneficial, as mesoscale air–sea interactions in the mid-latitudes are reported to demonstrate a significant impact on model representations of large-scale climate (Minobe et al. 2008; Kirtman et al. 2012; Ma et al. 2017). Further developments are needed to explore these exciting issues.

Data Availability Statement

The reforecast data used in this paper can be obtained from the Japan Meteorological Business Support Center, the Meteorological Research Consortium of JMA (the Meteorological Society of Japan), or the European Union's Copernicus Climate Change Service. The program code used in this study is not publicly available due to the management policy of JMA but may be available from the relevant authors for usage upon reasonable request, subject to permission from JMA.

Acknowledgments

CPS3 uses GSM and MRI.COM, which many staff members at JMA have been involved in developing. Besides the authors, T. Matsukawa, A. Nishimura, N. Otsuka, S. Ishizaki, and Y. Ichikawa also contributed to the development of CPS3. H. Kawai, M. Nakagawa, H. Yoshimura, K. Yoshida, Y. Takaya, and N. Oshima provided valuable advice on physical atmospheric processes, S. Urakawa consulted on SCUP, T. Toyoda advised on sea ice data assimilation, and T. Y. Tanaka provided advice on considering volcanic aerosols. Comments from two anonymous reviewers help improve the manuscript.

Footnotes

1A bug has been found that Typhoon Bogus was unintentionally excluded from JRA-3Q for the period after 2013; an updated version is being prepared at the time of writing. However, owing to the limited area and period affected by the bug, re-running of the forecast is not planned.

References
 

© The Author(s) 2023. This is an open access article published by the Meteorological Society of Japan under a Creative Commons Attribution 4.0 International (CC BY 4.0) license.
https://creativecommons.org/licenses/by/4.0/
feedback
Top