Back to Peter Oke’s HOMEPAGE




BODAS – Bluelink Ocean Data Assimilation System

Dr Peter R. Oke


Summary ● Description ● _Characteristics of BODAS


link to a complete list of publications


Oke, P. R., G. B. Brassington, D. A. Griffin and A. Schiller, 2007: The Bluelink Ocean Data Assimilation System (BODAS), Ocean Modelling, doi:10.1016/j.ocemod.2007.11.002.

Oke, P. R., A. Schiller, G. A. Griffin, G. B. Brassington 2005: Ensemble data assimilation for an eddy-resolving ocean model. Quarterly Journal of the Royal Meteorological Society, 131, 3301-3311.


The primary test-bed for BODAS in BRAN.






Alves, Oscar


Seasonal Prediction

Andreu-Burillo, Isabel


Coupled limited Are

Baird, Mark


Regional / Coastal ocean

Brassington, Gary


Global / Regional

Griffin, David



Gunson, Jim



Fiedler, Russel



Freeman, Justin



Huang, Xinmei



Hudson, Debbie



Jones, Emlyn



MacDonald, Helen



Mansbridge, Jim



O'Kane, Terry



Oke, Peter



Okely, Patricia



Schiller, Andreas



Sims, Holly



Sun, Chaojiao



Wedd, Robin



Yin, Yonghong






BODAS is an ensemble optimal interpolation system that uses an ensemble of intraseasonal anomalies from a free running model to estimate the background error covariances (BECs). The ensemble-based BECs are multivariate and inhomogeneous and are shown to reflect the length-scales, the anisotropy and the covariability of mesoscale oceanic processes.



Analyses of sea-level η, T, S, and horizontal currents (u,v), are computed by solving the analysis equations,





K=(ρring operatorP)H*(H(ρring operatorP)H*+R)-1,




Click to view the MathML source


is the state vector; superscripts a, b, o and * denote analysis, background, observed and matrix transpose, respectively (here the background field refers to a model-generated estimate of the ocean state at the analysis time; sometimes called the first guess); K is the gain matrix; H is an operator that interpolates from the model grid to observation locations; ρ is a correlation function used to localise the ensemble-based BECs in P; R is the observation error covariance matrix; and the open circles denote a Schur, or Hadamard, product (an element-by-element matrix multiplication). This formulation of the analysis equations is the same as that of Houtekamer and Mitchell (2001).

Estimates of the BECs in (2) are given by




where n is the ensemble size and A is an ensemble of model anomalies. Each anomaly field consists of all model variables included in (3). The calculation of the anomalies for an EnOI system is critical to the performance of the scheme. They should be computed in such a way that the scales of variability and features represented by the anomalies resemble the dominant errors of the model. For example, for application to OFAM, where we seek to correctly reproduce the mesoscale variability around Australia, we expect that the errors of an individual forecast will be dominated by the errors associated with mis-placement of eddies. Section 4 shows that this turns out to be true, although it is also true that errors at other scales exist. We hope to address this compromise in the future. Consequently, for application to OFAM, the ensemble of anomalies are generated by calculating intraseasonal anomalies derived from the spin-up run:


Click to view the MathML source


where α is a scalar that can tune the magnitude of the covariances for a particular application (for BRAN1.5, α=1); Click to view the MathML sourceis the ith intraseasonal anomaly, defined here as the 3-day average of the model state minus the seasonal cycle of the model spin-up run. Using this approach, anomaly fields in the ensemble look like typical eddy fields (in the high-resolution region). Because the ensemble in (5) is essentially a time series of anomalies, care is taken when constructing the ensemble to ensure that the time series is stationary, with zero mean and no trend. For both BRAN1.0 and BRAN1.5, BODAS uses an ensemble size of 72, with one anomaly from every month of the last 6 years of a 13-year model spin-up.

The BECs in (4) are localised in the horizontal around each observation in (2) using the localising correlation function ρ. Elements of ρ are defined by the quasi-Gaussian function of Gaspari and Cohn (1992), after Houtekamer and Mitchell (2001). Localisation has been shown to reduce the effects of sampling error for applications of an EnKF (e.g., Hamill et al., 2001) and EnOI (Oke et al., 2007). The localising correlation function in ρ forces the BECs to reduce to exactly zero, over L_ from an observation location. Note that in the application, we do not localise the covariances in the vertical direction. The present implementation of BODAS uses a uniform horizontal radial distance L = 8_, corresponding to an e-folding scale of about 3.5_. This has several ramifications for the system’s performance. Firstly, the rank of the estimated BECs in P is increased significantly. Using an ensemble size of n, the rank of P without localisation is at most n-1. By contrast, using an ensemble of only 72 with localisation with L = 8_ we estimate that the effective rank of ρring operatorP becomes Click to view the MathML source(Oke et al., 2005). This enables the assimilation system to determine analysis increments that fit the background innovations, given by (wo-Hwb) from (1), reasonably well. There are however, a few drawbacks to localisation. For example, analyses are not as dynamically balanced as they would be without localisation ([Mitchell et al., 2002] and [Oke et al., 2007]). Also, the inversion of the innovation covariance matrix (H(ρring operatorP)HT+R) becomes expensive, since the techniques for computational efficiency described by Evensen (2003) are not suitable when localisation is used. This makes the practical implementation of BODAS, described in the following sections, a technical challenge.

The primary test-bed for BODAS is BRAN. BRAN is a multi-year integration of OFAM, where BODAS is used to sequentially assimilate observations once every 7 days. Observations assimilated in BRAN include SLA from satellite altimetry, satellite-derived sea-surface temperature (SST) and in situ T and S profiles. Conceptually, BRAN is a three-dimensional time-varying synthesis of oceanic observations that uses OFAM as a dynamic interpolator. The sequential nature of the assimilation used here means that BRAN can be regarded as a series of 7-day “forecasts”. Of course, the skill of BRAN should exceed that of an equivalent real-time forecast system with a 7-day update cycle, because of the use of real-time surface fluxes and observations; and particularly because of the latency of real-time altimetry in the operational system. Despite these differences, we assess the short-range predictive skill of the system for SLA and SST in BRAN as a way of measuring the potential performance of the operational system.


Characteristics of BODAS

As stated above, one of the advantages of EnOI is that the BECs are inhomogeneous and anisotropic, reflecting the variability and length-scales of the ocean circulation. An example of the anisotropy and inhomogeneity of the ensemble-based BECs used by BODAS is shown in Fig. 4. The chosen examples show the localised ensemble-based cross-correlation between sea-level at a reference location and sea-level in the surrounding region. These fields demonstrate the region of influence of an observation at the reference location. Where the correlation is positive, the increment that is given by the term K(wo-Hwb) in (1), due to each observation has the same sign as the background innovation (wo-Hwb). So if the observed sea-level is higher than the background field at the reference location, in the absence of other observations in the same region, the solution to (1) will produce an increment that is also positive. The magnitude of the increment depends on the relative magnitudes of the estimated background error and observation error covariances (i.e., the relative magnitudes of H(ρring operatorP)HT and R). The structure of the increment also depends on the structure of the localised background error covariance, ρring operatorPHT.


Display Full Size version of this image (30K)

Fig. 4. Examples of the ensemble-based cross-correlations between sea-level at a reference location, denoted by the star, and sea-level in the surrounding region for a reference location (a) on the continental shelf, (b) over the continental slope and (c) over the deep ocean off eastern Australia (panel (d)). Contour intervals are 0.2; zero is bold, dotted is negative, correlations above 0.6 are shaded.

The reference locations for the examples presented in Fig. 4 are at 32.5_S, corresponding to the typical separation point of the EAC (Godfrey et al., 1980). These examples include correlations when the reference location is on the continental shelf at 115 m depth, where the correlation field has short decorrelation scales in the across-shore direction and long decorrelation scales in the alongshore direction (Fig. 4a). The long length-scales in the alongshore direction are probably a reflection of the covariability associated with northward propagation of slow-moving coastal trapped waves and the along-shore advection of the EAC and wind-driven circulation. Also shown is a correlation map when the reference location is over the continental slope at 1800 m depth (Fig. 4b), less than 50 km east of the shelf example (Fig. 4a). This field shows a more isotropic correlation, but with a tendency to have higher correlations to the east and south-east, in the typical direction of the EAC as it separates from the coast. The final example considered here shows correlations when the reference location is over the deep ocean at 4500 m depth (Fig. 4c). This field shows a somewhat isotropic correlation field, but also shows areas of weak negative correlation to the north-east and south-west. These negative correlations may be related to variations in sea-level associated with the typical eddy field in this region. The correlation fields presented in Fig. 4 highlight the anisotropy and inhomogeneity of the ensemble-based correlations, with relatively long quasi-isotropic correlations in the deep ocean, transitioning to very anisotropic correlations nearer the coast.

An example of the multivariate nature of BODAS is demonstrated in Fig. 5 and Fig. 6, showing the increments from a single observation analysis (i.e., where an analysis is computed by combining a single observation with a background field). In this example, an observation of sea-level at the coast is presumed to be 20 cm lower than the background sea-level. BODAS is used to calculate the increments to the full model state (η,T,S,u,v) in the surrounding region. We choose to use sea-level from a coastal location off South Australia at 140.35_E and 37.85_S. This location is arguably the best region for observing wind-driven, coastal upwelling along the Australian coastline (e.g., [Lewis, 1981] and [Kampf et al., 2004]). We therefore expect the increments to be consistent with a conceptual model of upwelling (i.e., consistent with the oceanic response to south-easterly winds).


Display Full Size version of this image (53K)

Fig. 5. Increments to sea-level (grey scale) and surface currents (vectors) based on a single observation of sea-level, denoted by the circle times operator, such that the observed sea-level is 20 cm lower than the background sea-level. The inset shows the region of interest off south Australia.


Display Full Size version of this image (89K)

Fig. 6. Increments to (a) sea-level, (b) across-shore currents (dashed contours are offshore; bold contour is zero) and (c) along-shore currents (solid contours are into the page; bold contour is zero); and the background field (solid contours) and analysis (dashed contours) for (d) T and (e) S, based on a single observation of sea-level at the coast as in Fig. 5.

Decreased sea-level at the coast may be due to a number of factors that include, but are not limited to, wind-driven upwelling. However, the dominance of locally wind-driven circulation in this region leads us to the expectation that it will also dominate the statistical properties of the ensemble fields here. We find that the increments to sea-level and surface currents in Fig. 5 are consistent with wind-driven upwelling, with reduced sea-level along the coast, a relatively strong coastal jet and weak south-westerly flow in the deep ocean. Fig. 6 shows the impact of the observation along a shore-normal section offshore of the observation location. This figure shows that the negative sea-level increment is strongest at the coast as we expect, and approaches zero moving offshore in a quasi-exponential fashion. Increments to the across-shore currents show an offshore flow of up to 5 cm s−1 over the top 30 m, that is consistent with an offshore wind-driven Ekman layer, and a weak, shoreward return flow at depth. Increments to the along-shore currents are consistent with a baroclinic, wind-driven coastal jet, with the strongest currents of 0.5 m s−1 at the surface near the coast. Fig. 6d and e shows the background and analysed T and S fields. Both of these fields show an uplift of isotherms and isohalines, with an implied vertical excursion of about 25 m near the shelf break and 40–50 m over the upper slope. From this analysis, we conclude that the statistical properties of the ensemble are consistent with a conceptual model of wind-driven upwelling in the region considered.

The analysis described above indicates that if the model does not produce an upwelling event, due to incorrect surface forcing for example, and we assimilate a single observation of sea-level at the coast that reflects the upwelling through reduced sea-level, then BODAS will compute increments that more closely match the observation at the coast, and in a manner that is consistent with wind-driven upwelling. While this feature is presented here as a benefit of EnOI, there is also a down-side. Suppose the model-observation mismatch is due to a mis-represented process other than wind-driven upwelling. Say, the encroachment of an eddy, or the propagation of a coastal trapped wave. Then the EnOI-derived increments will still be consistent with upwelling as in Fig. 5 and Fig. 6. This may not be desirable, and may result in a reanalysed state that is somewhat inconsistent with reality. However, we note that the example described here is very idealised. In practice, we typically assimilate multiple observations of various types (e.g., sea-level, T and S) and in the presence of additional observations that provide a more complete picture of the ocean state, BODAS will compute increments that more appropriately represent the true circulation.

Another limitation of EnOI that can be seen from the example described by Fig. 5 and Fig. 6 is that it implies a symmetry in the increments for an observation-model mis-match of the opposite sign. Suppose that the observation-model difference considered above is reversed. That is, the observation is 20 cm higher than the background sea-level. In this case, the increments simply have the opposite sign to those presented in Fig. 5 and Fig. 6. This might occur if a modelled upwelling is too strong, or perhaps if the model fails to represent a downwelling event. However, we note that studies into the dynamics of idealised wind-driven upwelling (e.g., Allen et al., 1995) and wind-driven downwelling (e.g., Allen and Newberger, 1996) on the continental shelf have shown that they are very different, and do not simply result in symmetric anomalies about some mean field as the EnOI-based increments would imply. These problems are the result of assuming ergodicity in the BECs. A better approach is to employ some flavour of the EnKF (e.g., [Evensen, 2003] and [Sakov and Oke, in press]), where the anomalies in (5) evolve in time and more accurately reflect the dominant dynamical processes for a particular time. However, these types of filters require an ensemble of states to be evolved, which is currently too computationally expensive for applications as large as OFAM.



[CMAR Home]

Last updated 22/09/06  | Legal Notice and Disclaimer | Copyright