Multi-Model Communication and Data Assimilation for Mitigating Model Error and Improving Forecasts

Yian CHEN Samuel N.STECHMANN

(Dedicated to Professor Andrew J.Majda on the occasion of his 70th birthday)

Abstract Models for weather and climate prediction are complex,and each model typically has at least a small number of phenomena that are poorly represented,such as perhaps the Madden-Julian Oscillation (MJO for short)or El Ni?o-Southern Oscillation (ENSO for short)or sea ice.Furthermore,it is often a very challenging task to modify and improve a complex model without creating new deficiencies.On the other hand,it is sometimes possible to design a low-dimensional model for a particular phenomenon,such as the MJO or ENSO,with significant skill,although the model may not represent the dynamics of the full weather-climate system.Here a strategy is proposed to mitigate these model errors by taking advantage of each model’s strengths.The strategy involves inter-model data assimilation,during a forecast simulation,whereby models can exchange information in order to obtain more faithful representations of the full weather-climate system.As an initial investigation,the method is examined here using a simplified scenario of linear models,involving a system of stochastic partial differential equations (SPDEs for short)as an imperfect tropical climate model and stochastic differential equations (SDEs for short)as a low-dimensional model for the MJO.It is shown that the MJO prediction skill of the imperfect climate model can be enhanced to equal the predictive skill of the lowdimensional model.Such an approach could provide a route to improving global model forecasts in a minimally invasive way,with modifications to the prediction system but without modifying the complex global physical model itself.

Keywords MJO,Multi-Model communication,Data assimilation,Kalman filter algorithm

1 Introduction

In this paper,in brief,we investigate an approach to weather and climate forecasting that uses multiple dynamical models.Such an approach becomes more relevant as prediction models have advanced from being local,regional models to being global models,since global models aim to represent the entire weather-climate system as a whole,but doing this accurately with a single model is difficult.

Furthermore,predictions are advancing beyond the short-term forecasts of only a few days,and now a major challenge is to make forecasts at the intersection of weather and climate,at lead times of weeks,months,or seasons in advance (see [5,21,39,50,54,69,78,82,94,95,97]).Such subseasonal to seasonal (S2S for short)predictions are related to some of the strongest signals of intraseasonal and interannual variability,such as monsoons,the Madden-Julian Oscillation (MJO for short),El Ni?o-Southern Oscillation (ENSO for short),sea ice,etc.,and these phenomena can be global themselves or linked globally through various teleconnections.

While global models such as global climate models (GCMs for short)are powerful tools for making predictions,they can have significant deficiencies,and the deficiencies cause prediction skills that are below the estimated predictability limits.For instance,GCMs commonly have deficient representations of the MJO (see [34,36,58,84]),and their predictions of the MJO are likewise below the estimated theoretical predictability (see [25,49,72,95,97,99]).Other phenomena that have deficient representations in GCMs include the intertropical convergence zone(ITCZ for short)(for which there is a serious double-ITCZ problem)(see[35,57])and the El Ni?o-Southern Oscillation (ENSO for short)(see [27,37]).

Improving a complex dynamical model such as a GCM is a very challenging task.Sometimes improving one aspect could cause a worsening of another aspect of the simulation.For example,such an issue has been noticed in attempts to improve the MJO in GCMs,where improvements in the MJO have sometimes been seen to be accompanied by a worsened climatological mean state (see [1,48]).Such attempts to modify GCMs to improve their simulations of the MJO are undoubtedly worthwhile.Nevertheless,the difficulties of such a challenging task can also be taken as motivation for alternative methods for simulating and predicting the MJO (as well as other challenging phenomena such as ENSO and sea ice).

Besides global models such as GCMs,an alternative approach to forecasting is to use simple models.For example,the MJO could be crudely represented by two indices,x1(t)and x2(t),which represent the amplitudes of the two empirical orthogonal functions (EOFs for short)of the MJO’s structure from principal component analysis (PCA for short)(see [102]).Then stochastic/statistical models can be used for predicting these types of indices (see [8,10]).Other types of simple models might incorporate aspects of atmospheric and/or oceanic fluid dynamics,thereby providing some added physical realism while still maintaining a simplified modeling approach with,say,O(100)degrees of freedom; examples of this type include the skeleton model for the MJO in [11,12,13,65,74,76,85,92,93] and models for ENSO in[91,108].Despite the simplicity of these types of models,they can often produce forecasts with significant skill.In fact,forecast skill with simple models sometimes meets or exceeds the skill of forecasts with comprehensive physical dynamical models such as GCMs (see [8,10,38,95,98,104]).Such results have been seen for not only the MJO but also other aspects of weather/climate variability,such as ENSO (see [80,81])and sea ice (see [7,59,100]).

Given the background described above,the main goal of the present paper is to investigate a multi-model communication strategy for forecasting with deficient models.As a concrete example,the case of the MJO will be used (although the general idea should be applicable to other scenarios described above as well).It is assumed that two models are available:(i)a deficient tropical climate model,which is meant to be an adequate model of climate but deficient in its representation of climate variability (the MJO),and (ii)a low-order model for the MJO,which has no representation at all of the rest of the weather-climate system.The goal is to utilize both models together to achieve the best possible forecast of tropical weather/climate variability,and to do so in a way that will emphasize each model’s strengths and mitigate each model’s deficiencies.

To carry out the multi-model communication strategy,the models will be allowed to assimilate data from each other while in forecast mode.In other words,each model will be run forward in time for a forecast,and occasionally each model will consider the other model’s forecast to be providing an“observation”,and this data will be assimilated.Schematically,the strategy could be described for two models as follows:

Here,u and v represent the state vectors of the two models,F andrepresent the dynamical operators,σ andrepresent(potentially random)forcing,the subscripts(0,1,2,···)represent points in time,the downward arrows represent evolution in time,and the horizontal arrows indicate the exchange of information via data assimilation.

For the specific scenario described above,the deficient tropical climate model will occasionally assimilate information about the MJO from the low-order MJO model forecast.It is hoped that the assimilated forecast data will help the deficient climate model to,e.g.,maintain and propagate an MJO with an enhanced skill that would otherwise be impossible if the climate model were run on its own.Such a strategy is plausible given that reanalysis products (see[40,e.g.])(which assimilate actual observational data into deficient climate models in order to provide a best estimate of the actual atmospheric state and its evolution)are commonly assumed to have reasonable (although not perfect)representations of the MJO and used to study its behavior (see [47,e.g.]),whereas free-running climate models that do not assimilate observational data are commonly seen to have deficient MJOs (see [34,36,58,84]).Here,it will be the low-order MJO model forecast that provides the “observations” of the MJO that could potentially drive a reasonable MJO in the deficient climate model’s forecast.On the other hand,in the opposite direction,one can imagine that the low-order MJO model would benefit from assimilating data from the climate model,which could provide information about the rest of the climate system and its influence on the MJO.

If such a strategy could be successful,it could allow one to leverage advancements in simplified models and predictions (see [8,9,10,65,74,76,92])to improve predictions with comprehensive global models.Furthermore,such an approach will potentially allow improvements in a minimally invasive way,since it is the prediction system that would be modified rather than the dynamics of the climate model itself.

Several simplifications will be utilized in the present paper.The main simplification will be the use of linear,Gaussian models for both the deficient climate model and the MJO model.On one hand,such a simplification is crude; on the other hand,it allows the use of the Kalman filter for the inter-model data assimilation,which will help provide a more transparent view of the workings of the methods here.Another simplification is that the deficient tropical climate model will use a crude vertical structure and no meridional structure(so it represents variability directly above the equator at y =0).On one hand,such a simplification will restrict the variety of modes of climate variability; on the other hand,this restriction can be taken as one way in which this climate model is deficient.

Another approach that involves the use of multiple models is the multi-model ensemble mean,which has been in use for some time (see [17,28,42,51,52,53,103]),but it involves multiple models only in post-processing,not through any type of direct communication between models.In allowing the coupling of multiple models,the multi-model communication approach is,to some extent,reminiscent of the topic of synchronization of chaotic dynamical systems(see [16,106])and the approach of superparameterization methods which couple the dynamics of macro-scale and micro-scale models (see [26,41]).Also,there has been some use of artificial intelligence,machine learning,neural networks,multifidelity methods(see[79]),and other data science methods for weather and climate applications(see[68,83]).To an extent,the method of the present paper could be viewed as a type of machine learning,where one model learns of some components of the weather-climate system from another model.While this idea is implemented here through a data assimilation algorithm,it is possible that other machine learning techniques could potentially be investigated for similar purposes in the future.

This article is dedicated to Andrew J.Majda to celebrate his 70th birthday,and it has drawn on his influence through his many contributions to topics of models for tropical weather and climate variability (see [2,4,44,45,61,65,105]),new convective parameterizations for GCMs (see [14,15,22,23,24,46]),data assimilation and filtering (see [6,9,19,20,29,30,31,32,55,62,64]),predictions of tropical weather and climate variability (see [8,10]),and multi-model ensemble strategies (see [3]),to name only a few.

The paper is organized as follows.The observational data and data analysis methods are described in Section 2.The two models-the deficient tropical climate model and the simple MJO model-are described in Section 3 and Section 4,respectively,including descriptions of the models themselves as well as their individual skill as models for data assimilation and forecasting.Finally,the model communication strategy is then investigated in Section 5,followed by conclusions in Section 6.Many details of the methods are presented in the appendices in order to streamline the presentation in the main text.

2 Data Source and Data Processing Methodology

Our tropical climate model (which will be presented below in Section 3)involves three variables:Wind (zonal velocity)u,potential temperature θ,and a variable q that represents moisture or water vapor or cloudiness.To assess the model’s skill,observation-based estimates of the three model variables are needed.This section provides a brief overview of the data sources and data processing methodology,and a more detailed description is given in Appendix 7.

Two data sources are used here.First,the National Oceanic and Atmospheric Administration (NOAA for short)interpolated outgoing longwave radiation (OLR for short)is used (see[56]).OLR is commonly used as a surrogate for cloudiness,and it is also somewhat well correlated with moisture anomalies (see [47,e.g.]).For the purposes of the present paper and its simple and deficient model,the distinction between moisture and cloudiness is not a major consideration,although in other frameworks,such as the MJO skeleton model,such a distinction is important (see [65,66]).Second,National Centers for Environmental Prediction-National Center for Atmospheric Research (NCEP-NCAR for short)reanalysis daily zonal winds and geopotential height are used (see [40]).Despite the fact that reanalysis data are a combined product of observational data and model dynamics,we sometimes refer to it as observational for simplicity.Both datasets have a horizontal spatial resolution of 2.5°×2.5°and a daily temporal resolution from 1 January 1979 to 31 December 2011.Since the connections between the model variables and observations are largely based on previous work (see [13,75,88,90]),further details are relegated to Appendix 7.

Here the influence of the seasonal cycle is removed via subtracting the time mean and first three harmonics of the annual cycle,based on the 1979-2011 observation record.And to further remove low-frequency variability,a 120-day-mean of the previous 120 days is subtracted,as is recommended by the Climate Variability and Predictability (see [96]).

Figure1 Plots of dimensionless (a)u,(b)θ,(c)q variables derived from observational data.The time period is 1 Jan 1979-31 Dec 1979.

An example of the observed evolution of the atmosphere is shown in Figure 1.Inspection of Figure 1 reveals some features occurring in MJO,especially the typical eastward-propagating waves,which are most visible in the wind u and also in the moisture/cloudiness q.

3 Deficient Tropical Climate Model

Our main goal of this paper is to illustrate an idea for multi-model communication,using inter-model data assimilation,as a strategy to improve forecasts.The improvements will come from mitigating the effects of model errors.As an illustration in a simplified setting,we utilize a simplified tropical climate model that has some significant deficiencies (Subsection 3.1).We intend to treat it as a complex climate model that can simulate the outlines of the climate but is deficient in its description of some important climate phenomena,such as the MJO.To finish the section,we describe the model’s free-running climate and variability(Subsection 3.2)as well as the model’s performance in data assimilation and prediction (Subsection 3.3 and Subsection 3.4,respectively).In later sections,we will use this deficient climate model as part of the multi-model communication strategy to improve its forecasting performance.

3.1 Model description

The deficient tropical climate model(DTCM for short)is a three-variable system of stochastic partial differential equations (SPDEs for short)with the following form:

where u is the zonal wind velocity anomaly,θ is the potential temperature,q is the water vapor anomaly andis spatiotemporal white noise.The equations have been nondimensionalized using the standard equatorial reference scales shown in Table 1 (see [60,88,89,90]).

Table1 Model parameters and reference scales.

In the model in (3.1)-(3.3),the left-hand side terms represent a linearization of hydrostatic atmospheric dynamics (see [89,e.g.]).On the right-hand side,represents momentum damping (Rayleigh friction),and -represents Newtonian cooling.The termrepresents a heat source and moisture sink associated with deep convection and precipitation.Finally,b?xxq and D*represent turbulent advection-diffusion of water vapor,represented as eddy diffusion plus stochastic forcing.

To put the model in perspective,it can be compared with other models that are similar in spirit.For instance,it is similar to convective adjustment models (see [18,70,71,87,107]),although a nonlinear switch is not used here(for simplicity),and eddy diffusion of moisture has been added here (for allowing a somewhat realistic representation of the background spectrum of tropical convection).From another perspective,one could view this as an extension of a model for stochastic advection-diffusion of tropical moisture(see[33]),where the present model has been extended to include equations for u and θ in order to represent waves.Only a single moisture variable has been included here,as opposed to a similar model that distinguishes between lower and middle tropospheric moisture (see [86]); and the use of a single moisture variable here will lead to deficiencies in the present model’s representation of the MJO and convectively coupled equatorial waves (CCEWs for short).

The values of the parameters are chosen to obtain the best(although quite deficient)match of the climatological statistics between observational data and free-running model simulations,as described in more detail below in Subsection 3.2.Such an approach for parameter estimation is somewhat common (see [33,86,e.g.]).The parameter values are summarized in Table 1.

The model only has two independent coordinates,x and t,which represent the longitude and the time,respectively.Therefore the model describes the tropical climate along the equator(y = 0)for all longitudes.The dependence on the vertical (z)coordinate would come from associating u,θ and q with vertical basis functions that represent the first baroclinic mode,as is common practice (see [65,88,89,e.g.]).

To facilitate comparison between the model and the observational data,the same spatial resolution should be used in the model.As described above in Section 2,the observational dataset has a spatial resolution of 2.5°×2.5°.Since the circumference of the Earth at the equator is 40132.8 km (26.7552 in nondimensional units),there are in total N = 144 data points along the equator,and the grid spacing is Δx=278.7 km (0.1858 in nondimensional units).

The model is a system of stochastic partial differential equations as functions of x and t,which are linear with constant coefficients and can be solved using Fourier transforms (see[33,62,63,86]).Some details of the solutions methods are described in Appendix 8.In the following paragraphs we refer to the model as the DTCM.

3.2 Free-running simulations

In this section,we present free-running simulations of the DTCM.The simulations can be compared versus observational data to assess the deficiencies in the model,and the comparison also provides a way to choose model parameter values (in order to obtain the best possible comparison of climatological statistics).

A free-running model simulation is shown in Figure 2.The plots are set up for comparison with the observational data in Figure 1.Both Figure 1 and Figure 2 represent only a small section from statistical steady states,and therefore the weather cannot be compared on a day-by-day basis,but the variability of the two datasets can be compared in a statistical sense.

Figure2 Plots of free-running model simulation of (a)u,(b)θ,(c)q variables.The time evolution is sampled one time every 8 hours.A one-year free-running simulation is plotted.

Note the deficiency in the waves of the climate model.In particular,note that DTCM appears to have equal amounts of waves that propagate in the eastward and westward directions(see Figure 2),whereas in observational data the eastward-propagating waves clearly dominate over the westward-propagating waves (see Figure 1).(One can identify eastward-propagating waves in these figures as streaks that are tilted from lower-left to upper-right in the plot.Conversely,one can identify westward-propagating waves in these figures as streaks that are tilted from lower-right to upper-left in the plot.)In addition to this visual appearance,this behavior can be quantified statistically by considering the power spectral density,as shown in Figure 3 and Figure 4.For observational data,the power spectra show larger power for positive wavenumbers(k ＞0),which corresponds with eastward-propagating waves(as the phase speed isand power is plotted for ω ＞0).In contrast,for the deficient climate model,the spectral power has equal magnitude for positive (k ＞0)versus negative (k ＜0)wavenumbers,and therefore has equal magnitude for eastward- and westward-propagating waves,respectively.In short,the climate model is deficient in the sense that it fails to properly simulate the different wave types (MJO,CCEWs,etc.)and their different propagation directions.(This deficiency is also evident from the model equations,which do not include the Coriolis term,which would introduce east-west asymmetry and plays a major role in determining the propagation of the MJO and CCEWs.)

Despite its deficiency in simulating waves,note that the model does at least present a good simulation of the climate variables in two ways.First,the model and observations have similar order of magnitudes of the anomalies (roughly 0.1 to 0.3,as seen by comparing Figure 1 and Figure 2).Second,the model and observations have the same broad structure of the power spectral density (see Figure 3 and Figure 4),which,to a first approximation,resembles the power spectrum of spatiotemporal red noise (see [33])(although the model is deficient in that it looks like the “background spectrum” and lacks the detailed additional shape imparted by the MJO and CCEWs (see [73,86,101])).

Figure3 (a)Power spectral density plot of observational q data.(b)Power spectral density plot of model-simulated q data.

The comparisons between model and observations (Figure 1 through Figure 4)were also used to determine the parameter values in the model.As a starting point,parameter values similar to [86] were used,and then slight modifications were made in order to obtain the best comparison between model and observations.Note that θ is noisier and more difficult to filter or predict accurately than the other two variables,so we mainly focused on fitting u and q better.This method of parameter tuning is somewhat ad hoc,although it is also quite common,and it serves the purposes of the present paper.

3.3 Filtering

In this section,we describe the Kalman filter method that is used for data assimilation,which will be used in later sections to provide the initial conditions for forecasts.The Kalman filter is an algorithm which uses two estimates together-an estimate from the physical model and an estimate from observational data-to produce a best estimate of the state of the system(u,θ,q).The details of the Kalman filter and the method of applying it to DTCM can be found in Appendix 8.

To assess the DTCM’s performance as a filter model,we use the following approach.We treat the observational data as a “true signal”,and then we add artificial observation noise to the “true signal” to get artificial “corrupted observations”.The Kalman filter is then applied using the artificial “corrupted observations” as one of the inputs.The strength (standard deviation)of the artificial observation noise is chosen to be 15% of the climatological standard deviation of the observational data (“true signal”),similar to other previous studies (see [8,e.g.]).

An example of filtering performance of q in Fourier wavenumber 1 is shown in Figure 5(a),and overall RMS error comparison for each longitude is shown in Figure 5(b).(Notice that the corrupted observations in panel (a)appear to be quite close to true signal,despite the presence of 15% noise in the corrupted observations; the closeness is due to the fact that this figure illustrates the signal for Fourier wavenumber 1,whereas the 15% noise is added at each grid point in physical space; so the impact of the noise on Fourier wavenumber 1 will appear smaller than the noise’s impact in physical space.)As a baseline for comparison in panel(b),the corrupted observations are seen to differ from the true signal with a RMS error of approximately 0.01.The filter result,on the other hand,appears to provide an estimate of the true signal that is actually slightly worse; its RMS error is roughly in the range of 0.010 to 0.013.Hence the DTCM with Kalman filter does not help us to get closer to the true signal,which is another indication of the model’s deficiencies.

Figure5 (a)Example of SPDE filtering compared with true signal and (corrupted)observations.The data we show in (a)is Fourier mode 1 of variable q with time period 1 Jan 1979-10 Apr 1979.(b)The averaged RMS error of (corrupted)observations and filter result for each longitude.The RMS error data are obtained by averaging independent experiments of time period 1 Jan 1979-24 May 1984.

3.4 Forecasting

In this subsection,the forecasting performance of DTCM is tested.The criterion we use to assess the prediction performance is the model’s predictability of two RMM indices.RMM 1 and 2 indices,whose full name is Real-time Multivariate MJO (RMM for short)series 1 and 2,are first introduced by [102].They are a pair of principal component(PC for short)time series based on projecting the observational data onto a pair of empirical orthogonal functions(EOFs for short).In many previous studies,the leading two EOFs generally appear as a pair and they are shown to be useful and skillful indices of monitoring the MJO.The spatial patterns of atmospheric variability,especially the eastward propagating signal,are well captured by RMM 1 and 2.Here we use these indices to aid in assessing the model’s ability of simulating and forecasting the MJO.

To derive RMM indices,EOF analysis is applied to the observational data.Here the combined field of (u,q)observational data is used to determine the EOF basis.This is essentially the same choice as in [102],except the step of normalizing each field by its global variance has been neglected here because,in our model,each variable has been nondimensionalized by its reference scale and therefore the nondimensional data for u and q already have very similar global variance.The seasonal cycle and further interannual variablity have been removed in data processing procedure.The spatial structure of two leading EOFs are shown in Figure 6.The variance explained by EOF 1 and 2 is 12.02%and 11.87%,respectively.Therefore,leading pair of EOFs are used to be low-dimensional approximations which desire to capture a considerable amount of variance.By projecting u,q combined field to leading EOF basis,we can get RMM 1 and 2 time series.

Figure6 Spatial EOF structure of leading two EOFs.

The forecasting method follows standard procedures.The Kalman filter described above is used for obtaining the initial conditions for the forecast.Ensembles of DTCM simulations can be run forward in time to provide the forecast.

The MJO prediction skill (pattern correlation)as a function of prediction lead time is shown in Figure 7.Each data point is averaged among independent experiments of time period 1 Jan 1979-24 May 1984.As we can observe in the figure,the prediction skill of DTCM decays dramatically with the lead time.As a baseline,we consider that if the averaged pattern correlation is greater than 0.5,we deem that the model can provide effective predictions under this given lead time.Therefore,DTCM can only provide effective predictions within lead time about 3 days.The prediction skill decays to nearly 0 after lead time 6 days.There is little difference in the prediction skill of the two indices RMM 1 and 2.

Figure7 Prediction skill of the deficient tropical climate model.Pattern correlation data are averaged among independent tests of time period 1 Jan 1979-24 May 1984.

4 A Simple Model for Intraseasonal Oscillations

4.1 Model description

With the stated DTCM as complex climate model which requires improvements,our next step is to construct a simple model which can capture the skeleton of the key features of MJO and use the simple model to improve the performance of DTCM.More specifically,since our main focus is on predicting RMM indices,the simple model should at least have better prediction performance for RMM indices.As stated in Subsection 3.4,RMM 1 and 2 indices are a pair of time series of leading EOFs which can express a considerable amount of atmospheric variability of MJO.Therefore,a complex-valued stochastic differential equation(SDE for short)is introduced in this section.The real part and imaginary part of the equation are aimed at representing the time evolution dynamics of the RMM 1 and 2 indices,respectively.In order to distinguish it from DTCM,we refer to the SDE model as MJO model in the following text.The complex-valued SDE has form

where γ,σ ＞0 and ω are real,dW(t)is a complex Gaussian white noise with form

Each component is an independent Gaussian white noise.The remaining problem is how to determine the parameters by the given observation time series.The idea is to choose parameters to let the MJO model simulation and the observational data have the same(or close to the same)statistics,especially the autocorrelation function.The MJO model is exactly solvable.Thus a formula of autocorrelation function in terms of the 3 unknown parameters can be derived.We can compare the formula with the climatological statistics to make them match as closely as possible.Detailed calculations can be found in Appendix 9.Here we choose γ =8.8,ω =12.1,σ =1.84,in nondimensional units; in dimensional units,these parameters are γ =0.088 day-1,ω =0.121 day-1,σ =0.184.Figure 8 shows that there is a reasonably good fit,for such a simple model,between the autocorrelation function of the observational RMM series and the SDE model simulation.

Figure8 Real parts of autocorrelation function of observational RMM indices (solid)and SDE simulated indices (dashed).

4.2 Free-running simulation

In this section,the connection is drawn between (i)the MJO model for the RMM indices and (ii)DTCM.The connection is made by utilizing the EOF structures to associate spatial structures u(x,t)and q(x,t)to the MJO model time series,and this connection will be used in the multi-model communication strategy in Section 5.

A free-running simulation of the MJO model will be used to illustrate the model’s variability and the reconstruction of u(x,t)and q(x,t).Let

Divide the equation into its real part and imaginary part,the equation (4.1)is also equivalent to the two dimensional system

The system is then solved numerically using the standard Euler-Maruyama method with a time step of 0.01 in nondimensional units.

Given time series of X1(t)and X2(t),we can reconstruct a corresponding combined field(u(x,t),q(x,t))by the determined EOF basis functions,which were presented earlier in Figure 6.An example of reconstructed field is shown in Figure 9.The result is similar to the OLR reconstruction shown in [102].Even though it is a simple model,the typical eastward propagation signal is well simulated.Through these reconstructed signals of u(x,t)and q(x,t),the MJO model can be connected with DTCM in the multi-model communication strategy in Section 5.

Figure9 Signals of (a)u and (b)q reconstructed from the MJO model simulated RMM 1＆2 time series,using the model in(4.3)-(4.4).The duration of the plotted data is 1 year.The duration of the plotted data is 1 year.

4.3 Filtering

Data assimilation will be used to provide the initial conditions for the MJO model forecasts(described below in Subsection 4.4).As a data assimilation strategy,since the MJO model is linear,we use the Kalman filter to estimate the values of the RMM indices.Before proceeding with forecasts,we assess the filter’s performance in the present section.

Figure10 Example of filtering comparing with true signal and (corrupted)observations.The data shown are RMM 1 index with time period 1 Jan 1980-19 Jul 1980.

The detailed filter setup can be found in Appendix 10.An example of filter performance is shown in Figure 10.The filter typically provides a smoothed signal,and it often misses the maximum amplitudes of the strongest events (such as the strong event in Figure 10 from day 60 to day 120).While the filter does not appear to offer an improved estimate of the state of the RMM indices,it will still be used to provide initial conditions for forecasts,since the idea is to create a setup that is similar to an actual operational forecast scenario.

4.4 Forecasting

The simple MJO model has much greater forecast skill than the DTCM,as illustrated in Figure 11.As we can observe in the figure,the MJO model can provide effective predictions up to lead time 5-6 days,which is much better than DTCM.Therefore,it is hoped that,in the model communication strategy presented below,the information from this simple MJO model will be able help improve the DTCM’s MJO forecasts.

Figure11 MJO model prediction skill of RMM indices 1 and 2.The pattern correlation data are averaged among 5000 independent tests.The testing time window is 1 Jan 1980-9 Sept 1993.

Note that the skill of this simple MJO model is less than other simple MJO models that have been presented in other studies (see [8,10,95]).Nevertheless,it is useful in that it is easy to formulate and simulate,and it will still serve the purposes of the present paper,since its forecast skill is greater than the DTCM’s forecast skill.

5 Model Communication for Improved Forecasts

In previous sections,two different types of models were described with different sets of pros and cons:The DTCM has certain defects and it has very limited prediction performance of the RMM indices,whereas the simple MJO model can only simulate two RMM time series but its prediction performance is much better.

Now,in this section,the key idea is to allow the two models to ‘communicate’ with each other during the forecast.By letting the comprehensive model assimilate information from the simple model,we want to improve the prediction performance of the comprehensive model to reach or even possibly exceed the performance of the simple model.More specifically,we want the DTCM to assimilate the good MJO predictions given by the MJO model.

Three main advantages of this method are as follows.The first one is that by assimilating data from simple model,we find a way to improve the comprehensive model and avoid to make any changes of the model,as comprehensive models may be too complex to modify and if modifying some parts of the model,it’s hard not to touch anything else.The second one is we now have more options to make specific improvement of our model.By choosing several targeted indices and constructing focused simple models,we are able to work on improve the model’s performance of these specific indices.The third one is we can also have controls of the data assimilating process.Choose observation errors to represent the degree of trust to the simple model and choose proper assimilation frequencies to involve more information or reduce computing costs.More details are discussed in the following subsections.

In this first investigation of model communication,a simple Kalman filter method is used to achieve communications between models.When making forecasts,we run the MJO model predictions of RMM indices independently,and then these RMM index predictions are considered“observations”by the DTCM.The DTCM assimilates the RMM index predictions to hopefully get better prediction results.The details about how to apply the filter algorithm in the process is provided in Appendix 12.

5.1 Forecast skill with model communication

The forecast skill with model communication is shown in Figure 12.Four cases of model communication are shown,and they are labeled by the strength of artificial “observational noise” that is added to “observed” RMM indices supplied by the simple MJO model.The four noise strengths are 0%,5%,10% and 15% of the climatological standard deviation of the RMM time series.It is reasonable to include some noise,partly because the MJO model’s predictions are themselves not perfect,and partly as a way to assess the impact of errors in the prediction system.The noise strength could also be taken as an indicator of the trustworthiness of the simple MJO model predictions,where the extreme case of 0% indicates that the MJO model predictions are trusted completely relative to the DTCM predictions of the RMM indices.While the noise strength is simply assigned a fixed value in the tests shown here,it would be interesting in the future to allow the different models to learn the trustworthiness of the other models’information based on the skill of past forecasts.Also shown for comparison in Figure 12 are the forecast skills of the MJO model (SDE model)and the DTCM (SPDE model).

Figure12 Prediction skills of (a)RMM 1 and(b)RMM 2.The options provided are MJO predictions (SDE),DTCM predictions (SPDE),and model communication performance with observation error 0%,5%,10%,15% of climatological standard deviation.

In brief,one can see in Figure 12 that the skill of model communication can become equal to that of the MJO model (SDE).Compared to the skill of the DTCM (SPDE model),this represents a significant improvement.While it is perhaps not surprising that the model communication skill is equal to the MJO model (SDE)skill in the case of 0% noise (since the forecast target is the same as the assimilated RMM index data,and the models here are both linear),it is interesting to note the changes in forecast skill for different noise strengths.In particular,relative to 0% noise,the case of 5% noise has essentially the same skill,and the case of 5%noise has only a somewhat small reduction.Furthermore,beyond these points of note for RMM 1,notice that RMM 2 is forecast skillfully via model communication for all considered noise strengths.Differences in forecast skill for RMM 1 versus RMM 2 have also been noted in some past studies with other models (see [95,e.g.]).Overall,these results suggest the possibility of increasing the DTCM’s forecast skill through a model communication strategy.

5.2 Choice of assimilation frequencies

In the results described above,each model communication result is obtained by assimilating“observations” every 1 day (in nondimensional units,Δt=3).Since the MJO model provides better predictions of the RMM indices,we would expect that assimilating information from the MJO model more often would generate better results.Also note that assimilation frequencies lead to no difference when the artificial observational noise is 0.In Figure 13,the observation error used is fixed to be 15%.Take RMM index 1 as an example.The pattern correlations of model communication results of different assimilation frequencies are shown in Figure 13.As we can observe in the figure,assimilation frequencies do affect the forecast skill of model communication.Higher“observation”/communication frequencies tend to provide better results but more calculation may be involved in the process.

Figure13 Model communication prediction skill of RMM 1 for different assimilation frequencies,all with relatively strong artificial observational error of 15%.

5.3 Comparison with multi-model mean

Another way to utilize multiple models in making a forecast is called the multi-model mean(MMM for short).To create the multi-model mean,each forecast model is run independently,with no communication between the different models.For instance,here the result would be a value xsimplepredicted by the simple MJO model alone,and also a value xDTCMpredicted by the DTCM alone.The multi-model mean prediction,xMMM,is then defined simply as the average:

The MMM prediction skill is shown in Figure 14,in comparison with the prediction skills of the DTCM (SPDE),the MJO model (SDE),and model communication (with 5% standard deviation artificial observation error).Multi-model mean prediction skill falls in between the MJO model and DTCM.It does provide some improvement over the DTCM by itself,and it does not require any communication between forecast models during the forecast stage.However the data assimilation-based model communication strategy we stated in Subsection 5.1 provides better results.

Figure14 Multi-model mean prediction skill (crosses)in comparison with the MJO model(SDE,solid),the DTCM (SPDE,dashed),and the DTCM with 5% observation error.

6 Conclusions

A model communication strategy was investigated for allowing multiple forecast models to assimilate information from each other in order to overcome individual model deficiencies.The results show that a deficient climate model’s skill at forecasting the MJO can be significantly enhanced through communication with a skillful,simple MJO model.The multi-model communication strategy is also able to outperform the multi-model mean,which utilizes multiple models only in post-processing and not in any type of direct communication during the forecasts.

In the future,it would be interesting to test the model communication strategy in a setup that is closer to an operational forecast.In the present study,a simplified setup was used in order to allow transparent evaluation of different components of the strategy.A future next step could utilize an actual global forecast model rather than the simplified linear Gaussian model used here,and it could utilize more sophisticated data assimilation methods for the intermodel communication rather than the simple Kalman filter used here.Also,the nonlinearity of an operational forecast model would cause the assimilation of MJO EOF information to be spread to other components of the weather-climate system other than the MJO,and it would be interesting to see whether such spreading of information is beneficial or detrimental to the forecast of other components of the weather-climate system other than the MJO.

7 Appendix:Data Processing

7.1 Estimate the first baroclinic component of u and θ

The zonal wind velocity data and geopotential height data are used to estimate u and θ,respectively.Two pressure levels,850hPa and 200hPa,are used to represent the bottom and the top of the troposphere,respectively.Then associate the data with contributions from two vertical modes,a barotropic contribution and a first baroclinic mode contribution.Take zonal wind velocity u as an example.The velocity of different pressure level can be expressed as

where uBTrepresents the barotropic contribution and uBCrepresents the first baroclinic mode contribution.Given the data from two pressure levels,the first baroclinic component can be estimated by

Similarly,the first baroclinic component of the geopotential height can be estimated as following:

And consider the hydrostatic balance,which has formIt results in the expression of the potential temperature anomaly θ from the model in terms of the geopotential height anomaly Z from the data

Since we use the standard reference scale to nondimensionalize the model,the reference scale for u and θ can be directly derived.The reference scales for u and θ are 52.08 m/s and 15.6 K,respectively.

7.2 Estimate of the water vapor anomaly q

The OLR data are used to estimate water vapor anomaly q.Since water vapor can absorb certain wavelengths of OLR adding heat to the atmosphere.High concentration of water vapor tends to reduce OLR below clear sky values.Thus we can use OLR data as simple indicators to estimate the water vapor anomaly

The reference scale for water vapor anomaly q is deduced by the following steps.In a recent study (see [88]),the water vapor reference scale is 6.27 g·kg-1.We use the same reference scale here.However unit conversions are needed to apply the reference scale to the OLR dataset above.The observational OLR dataset has unit wm-2.Firstly,we use a OLR-to-heating-rate conversion factor which is 0.06 Kday-1(wm-2)-1.This conversion factor is the same as in [88]which is based on a recent study in [90].Now the unit becomes Kday-1.We want to further convert it to mm·day-1.Assume an area of 1 m2,now we estimate the mass of air in the atmospheric column.Consider the pressure of air is decaying exponentially with the height and the scale height is Hp=10 km.The function of pressure in terms of height z can be expressed by

Note that the scale height depends on the air temperature,by turning the scale height to a function of air temperature we can get a better estimate of the pressure.But here we only want to get a rough estimate so a constant value of scale height is used in the computation.According to the ideal gas law,assume the temperature is constant,the density of air is proportion to the air pressure.Thus the density of air can be expressed by

Use ρ0= 1 kg·m-3.Integrate the density of air from z = 0 to infinity and times the area.An estimation of the mass of air in the column is 10000 kg.Assume a depth of 1 mm of liquid water (after the water vapor is condensed to liquid),we use that the density of water is 103kg·m-3and the latent heat factor of water is 2500 kJ·kg-1.If this amount of water is converted from water vapor to liquid,the heat release in the process is 2.5×106J.If this amount of heat is used to warm the air in the atmospheric column and we use the specific heat of air is 1004 J·(kg·K)-1,the temperature of air will raise 0.249 K.Thus we get the unit conversion relation:1 K·day-1≈4 mm·day-1.

Now the unit becomes mm·day-1.Since the reference scale we use is 6.27 g·kg-1and the 1 m2column has 10000 kg air,that equals to 62.7 kg water vapor,which also equals to 62.7 mm liquid water in the column.This value can only apply to the air in low altitude or in other words,near the surface.If taking the vertical structure of water vapor into account,a more accurate expression of the ratio is needed.By many other studies (see [43,71,77]),the water vapor ratio of the vertical atmospheric column also decays exponentially with the height.Thus the water vapor ratio can be expressed by

Here we also use Hq=10 km as scale height and q0=62.7 mm.Note that the reference scale is the ratio of water vapor to the total humid air.And the total amount of humid air decays with height according to(7.7).Thus the reference scale used in our model can be calculated by

To associate the water vapor reference scale with the observational OLR data,divide the reference scale by τc=96 h=3 days.Thus the heating-rate reference scale is 10.45 mm·day-1=2.61 K·day-1=43.54 wm-2.By using this reference scale,we can nondimensionalize the OLR data.

7.3 Meridional projections

By previous step,we find proper surrogates for all three variables(u,θ,q)in our model.Next step is to use the meridional mode decomposition method to capture the skeleton of the data.By utilizing the meridional basis functions which take forms of the parabolic cylinder functions.Each variable can be expressed as linear combinations of parabolic cylinder functions.Take uBCas an example,

where φm(y)are basis functions and umrepresents the m-th meridional mode.Use the 0 meridional mode u0(x,t)as the skeleton of the whole data/ This quantity can be obtained by the projection

This meridional projection also reduces the 2D (x,y)dataset to a 1D (x)dataset.Same formula can also apply to θ and q.We can project θBCand q to φ0(y)and get θ0(x,t),q0(x,t)respectively.We use (u0,θ0,q0)as an estimate of variables (u,θ,q)in our model.To sum up,we apply the above data processing method to OLR,zonal wind velocity and geopotential height dataset.Then nondimensionalize the data by their corresponding reference scales and project them to meridional basis functions to get the observation data we use in our model.The observation data of three variables derived from the dataset are shown in Figure 1.Eastward propagating zonal wind can be clearly observed from Figure 1(a).

8 Appendix:Solving the Tropical Climate Model

Define the discrete Fourier transform for u as following:

Analogous formulas can also be apply to θ and q.Then we get the expression of the system in Fourier space,

Now the system has been transformed to a three-dimensional linear system.Denote(uk,θk,qk)Tby Vk.The system can be rewrote by matrices form

Note that the coefficient matrix depends on the Fourier wavenumber k.But each Fourier mode is independent with others which makes it possible to focus on solving one system and apply the same method to all other Fourier modes.Next,diagonalization is performed to the system.Let Akbe the eigenvector matrix (whose column are eigenvectors)of Fkand Λkbe the eigenvalue matrix (whose diagonal entries are eigenvalues).And let Uk= A-k1Vk=(Uk,1,Uk,2,Uk,3)T.Then Uksatisfies the diagonalized system

Each variable can be solved analytically by integrating exactly and applying the properties of.The exact solution is

Note that the integral in (8.9)is a Gaussian random variable at each time t.The mean and variance can be easily derived.Then Uk,j(t)has mean eΛk(j,j)tUk,j(0)and variance

This method is called semi-analytical solution method because the solution is derived analytically which means no numerical integration error occur,but the solution contains random term which may cause sampling error in realization.This formula can be used to numerically realize the model.Each Fourier mode can be evolved forward independently.For Fourier mode k,in each time step,given the initial data vector Vk,n= (un,θn,qn)T,use the matrix Akto diagonalize the system and get Uk.Then draw a Gaussian random variable with mean and variance given in (8.10).Then use matrix Akagain to change back and get numerical solution Vk,n+1.Treat it as initial condition in next time step to move forward.

9 Appendix:Kalman Filter for the Tropical Climate Model

A detailed description of the algorithm is introduced in [67].In the process,two kinds of information is involved-an estimation from the previous time step and an estimation from current measurement.Assume u is an M-dimensional unknown variable which requires estimation.And assume the filter model has form

where F is a complex dynamic matrix and σm+1is an M-dimensional complex Gaussian noise vector with zero mean and covarianceThus we assume that the equation(9.1)can give an unbiased estimation of the current state um+1from the previous state.Two states of the unknown variable are provided in the process.The prior state um+1|mis an estimation of the true state prior to the knowledge of the observation at time m+1.The posterior state um|mis an estimation of the true state after considering the observation at time m.The relation is given by

And assume noisy observations satisfy

where G is the complex observation matrix and σ0mis an M-dimensional Gaussian noise vector with zero mean and covarianceKalman filter is a two-step algorithm.In the prediction step,estimate prior mean state and prior error covariance by

In the correction step,the kalman gain matrix is defined as

And the estimations of the posterior state are

The posterior state is used as the kalman filter estimate of the true state.To associate the filter algorithm with DTCM,several choices should be made.Note that the model can be linearized in Fourier space.The true signal we want to estimate is Vk= (uk,θk,qk)T,for each Fourier wavenumber k.Since there is an easy transition between Vkand Ukstated in Appendix 8,we will estimate Ukinstead.The time evolution dynamics of Ukare given by equation (8.9).For each time step,the dynamic matrix

And the model error covariance matrix

where Var[Uk,1(Δt)U*k,1(Δt)] is given by equation(8.10).And the observations we use are data derived from dataset in Section 2.Thus the observation matrix is

Note that the observation error covariance matrix is a measurement of the extent of how we trust the observation data.We choose 15% of the standard deviation of each variable from the observation data as observation error,as is common in other similar studies (see [9,32]).Thus the observation error matrix is

Apply Kalman filter algorithm to each time step and provide the filter results given by DTCM.

10 Appendix:Kalman Filter for Simple MJO Model

Implement the standard Kalman filter to simple MJO model.For consistency,we use the same notation as in Appendix 8.As the filter model,we use the complex MJO model (4.1),and we create complex-valued observational RMM indices by using

To choose the filter parameters such as the model’s evolution operator and noise magnitude,consider the analytic solution of the MJO model:

Such a formula will be used to evolve the model from one observation time to the next,with observation time interval Δt.Therefore,if using XTas the signal we want to estimate,the filter model dynamic matrix will be

The model error matrix

and the observation matrix GMJO= 1.To create a scenario with “observational error”,we treat the observed RMM indices as the “true signal” and add artificial observation error to get artificial“observations”.For this purpose,we use 15%of the RMM indices’ standard deviation as the “observation error” (R0=15%(std(RMM 1)+i×std(RMM 2))).

11 Appendix:Determining Parameters in the MJO Model

Autocorrelation function is a useful tool to analyze the repeating patterns in a signal.In order to determine the parameter values,we assume the RMM indices satisfy the form of the MJO model(4.1).We first derive the autocorrelation function of the given SDE equation.Note that the MJO model is a complex Ornstein-Uhlenbeck process.The analytic solution is

That means at any time t,the solution is a random Gaussian variable with mean and variance

Define the autocorrelation function

Then we can calculate the autocorrelation function of the SDE,

We denote the integral of the autocorrelation function as

Thus the integral of the autocorrelation function can be used to calculate the parameters.And furthermore,by equation (11.3),σ can also be calculated,

12 Appendix:Kalman Filter for Model Communication

To implement standard Kalman filter for model communication,we now use RMM time series simulated by the simple SDE model.The following equations summarize the relations between signals,

The true signal we want to estimate is (uk,θk,qk)and the “observations” are RMM time series (PC1,PC2).These two signals can be evolved by equations (12.1).However the true signal is not fully observed since only two components (uk,qk)contribute to the time series.In this case,the relation between true signal and observation can be expressed in equation (12.2).

Now the problem is to find observation matrix G and properties of the error term σ0m.To determine the observation matrix G,an alternative representation of time series in terms of Fourier modes is required.Firstly,by EOF analysis,the(u,q)combined field can be decomposed by the EOF basis.Let φ1,φ2,··· ,φNdenote the EOF basis vector and PC1(t),PC2(t),··· ,PCN(t)denote their corresponding time series.Note that PC1(t),PC2(t)are exactly RMM 1 and RMM 2 indices.Let u and q be of form time×longitude.Then we can write

Since EOF basis vectors are pairwise orthogonal,the PC time series can be calculated by inner products

Let u = (u1,u2,··· ,uN),q = (q1,q2,··· ,qN)and φj= (φj,1,φj,2,··· ,φj,2N)T.By discrete Fourier transformation,

Then plug in equation (12.4),

Therefore,by equation (12.7),we associate the EOF time series with Fourier modes.Consequently,we can bring the Fourier representation of the time series in the filter algorithm.According to the diagonalization in Appendix 8,Vk=AkUk.Denote

Then

Plug in equation (12.7),

With the above preparations,now we turn to the filter algorithm.For each time step tm,the true signal which requires estimation is um= (U0,1,U0,2,U0,3,··· ,UN-1,1,UN-1,2,UN-1,3)T.Note that we omit the time subscript m in the vector to avoid redundant subscripts.Note that the observations (RMM 1,RMM 2)are time series vm= (PC1(tm),PC2(tm)).According to equation (12.11),the observation matrix is

And assume RMM time series 1 and series 2 are decorrelated with each other.Then the observation error matrix is a diagonal matrix.Different choices of observation error are discussed in Subsection 5.1.