- Email: [email protected]
Assimilation of altimeter data into a quasigeostrophic ocean model using optimal interpolation and eofs
JOURNAL
ELSEVIER
OF
MARINE SYSTEMS Journal of Marine Systems 6 (1995) 125-143
Assimilation of altimeter data into a quasigeostrophic model using optimal interpolation and eofs
ocean
M.M. Rienecker, D. Adamec Code 971, Laboratory for Hydrospheric Processes, NASA /Goddard
Space Flight Center, Greenbelt, MD 20771, USA
Received 20 September, 1993; revised and accepted 24 March, 1994
Abstract An ensemble of fraternal-twin experiments is used to assess the utility of optimal interpolation and model-based vertical empirical orthogonal functions (eofs) of streamfunction variability to assimilate satellite altimeter data into ocean models. Simulated altimeter data are assimilated into a basin-wide 3-layer quasi-geostrophic model with a horizontal grid spacing of 15 km. The effects of bottom topography are included and the model is forced by a wind stress curl distribution which is constant in time. The simulated data are extracted, along altimeter tracks with spatial and temporal characteristics of Geosat, from a reference model ocean with a slightly different climatology from that generated by the model used for assimilation. The use of vertical eofs determined from the model-generated streamfunction variability is shown to be effective in aiding the model’s dynamical extrapolation of the surface information throughout the rest of the water column. After a single repeat cycle (17 days), the analysis errors are reduced markedly from the initial level, by 52% in the surface layer, 41% in the second layer and 11% in the bottom layer. The largest differences between the assimilation analysis and the reference ocean are found in the nonlinear regime of the mid-latitude jet in all layers. After 100 days of assimilation, the error in the upper two layers has been reduced by over 50% and that in the bottom layer by 38%. The essence of the method is that the eofs capture the statistics of the dynamical balances in the model and ensure that this balance is not inappropriately disturbed during the assimilation process. This statistical balance includes any potential vorticity homogeneity which may be associated with the eddy stirring by mid-latitude surface jets.
1. Introduction The assimilation of oceanographic data into numerical models is essential not only for predictive purposes, but also for diagnostic purposes, to help further the understanding of the ocean circulation on regional, basin and global scales. Given the sparsity and asynopticity of hydrographic and current meter data, oceanographers look forward to the availability of more numerous data through remote sensing. Satellites provide many measurements of several ocean variables;
but the measurement of the ocean surface topography by satellite altimeters provides the most direct connection to ocean dynamics since the altimeter signal is essentially a measure of vertically integrated mass. The problem arises as to how best to use this surface measurement to gain information about subsurface ocean variability. This paper describes a technique for projecting the surface information provided by the satellite altimeter to the ocean subsurface in an efficient manner and compares the results with a similar technique which relies entirely on the dynamical
0924-7963/9.5/$09.50 0 1995 Elsevier Science B.V. All rights reserved SSDI 0924-7963(94)00013-2
126
MM. Rienecker, D. Adamec /Journal
model to derive the subsurface variability consistent with that at the surface. Altimeter data has been incorporated into dynamical models by either assimilating the information only at the surface of the model (e.g., Hurlburt, 1986; Holland, 1989; Holland and Malanotte-Rizzoli, 1989; Verron, 1990, 1992; White et al, 1990a, b; Marshall and Berry, 1991) or by projecting the information in some manner in the vertical (e.g., De Mey and Robinson, 1987; Haines, 1991; Mellor and Ezer, 1991). The projection techniques have used either statistical information or dynamical “constraints”. The former methods include the use of empirical orthogonal functions of density perturbations determined from observations (De Mey and Robinson, 1987) or the use of regression techniques to estimate subsurface patterns (Mellor and Ezer, 1991). Holland and Malanotte-Rizzoli (1989) rely on model dynamics to project the surface information throughout the water column, viz., they assimilate along-track potential vorticity (i.e., the contribution to the potential vorticity from along-track derivatives of surface height) in the surface layer of their quasi-geostrophic (QG) model. Since the QG layer streamfunction is coupled vertically through the definition of QG potential vorticity, the surface information is passed immediately through the water column. In contrast, Haines (1991) inserts the surface streamfunction in the topmost layer only of the QG model and applies constraints on the potential vorticity in the subsurface layers and finds that the surface information is effectively projected through most of the layers. Different techniques are also used to incorporate the horizontal distribution of information. The most common methods are direct insertion, nudging and optimal interpolation. For direct insertion, the model streamfunction is simply replaced by the observation, as, e.g. in Haines (1991) who assumes that observations are available at all model grid points. For nudging an extra forcing term, proportional to the difference between the model variable and the observation, is added to the dynamical equations (e.g. Holland and Malanotte-Rizzoli, 1989). Optimal Interpolation (011, the data assimilation technique most com-
of Marine Systems 6 (1995) 125-143
monly used in operational numerical weather prediction (see, e.g. Bengtsson et al., 19811, produces an analysis field which is a linear combination of model-generated streamfunction and observed streamfunction. By using weights based on estimated error covariances for the model forecasts and the measurements the linear combination minimizes, in a least-squares sense, the errors in the analysis. White et al. (1990a, b), in identical-twin experiments, assimilated simulated Geosat altimeter observations into a 3-layer eddy-resolving, flat bottom, square basin QG model using 01 in the surface layer of the model only. Their first study assimilated altimeter differences into a model with initial fields 1% in error from the reference fields and found that the nonlinear regime of the mid-latitude jet limited the ability of the observations to constrain the model, so that the model had to be re-initialized periodically. Their second study placed stronger controls on the influence of the observations through the error covariance and in repeated use of the observations. The simulated altimeter observations were referenced to the model mean streamfunction and the study found that the assimilation procedure was effective in increasing the correlation between analyses and the reference ocean for all layers, however, the integrations were conducted for only 50 days. Marshall and Berry (1991) note that whereas surface information appears to be useful in constraining the subsurface flow in a two-layer model (depending on the satellite repeat time), this ability is severely limited in models with higher vertical resolution. It is apparent that techniques should be sought to the aid the model in its dynamical extrapolation. In this paper, an ensemble of fraternal-twin experiments is used to assess the assimilation, by optimal interpolation, of simulated surface altimeter data into a 34ayer square-basin QG model with bottom topography, forced by a wind stress curl distribution which is constant in time. The simulated data, along altimeter tracks with spatial and temporal characteristics of Geosat, are extracted from a time-series of surface-layer streamfunction in a reference ocean generated by
M.M. Rienecker, D. Adamec /Journal
the QG model with slightly different model parameters. The results from assimilation of the projection of the data onto empirical orthogonal functions (eofs) of the model streamfunction is compared with those from assimilation of the data into the surface layer only of the model. The use of model-generated eofs can be considered a compromise between the use of eofs derived by observations as in De Mey and Robinson (1987) and Dombrowski and De Mey (1992) and the use of a specific dynamical constraint, as by Haines (1991). The resulting analysis fields are therefore filtered according to the model’s dynamics and statistics. This paper is organized as follows: the assimilation system and the fraternal-twin experiment design is described in Section 2, the results are presented in Section 3 and conclusions summarized in Section 4.
2. The assimilation
system
2.1. The QG model The dynamical model used in the assimilation system is a three-layer QG model with solution technique as described by Cummins and Mysak (1988). For an N-layer model, the QG streamfunction is found as the solution to
+
fll j-+-l,2
-
wk+1/2
)
-44%
k
curl7 - $@v2*,
+ a,,--
7
PI+,
k=l,...,N;
(I)
and ;cs,
-
&+I)
=J(vh
-
h+*7
4&+1,2)
&+1/2 -
k=l,...,N-1.
---wk+l/2,
fo
(2)
of Marine Systems 6 (1995) 125-143
121
In the above, V2 is the horizontal Laplacian operator, J is the Jacobian operator, f = f0 + fly is the Coriolis parameter on the P-plane centered at 35”N, p is the density; g is the acceleration due to gravity; h, is the bottom topography; Hk are the average layer thicknesses, I?,, = 1, if 1= k, 0 otherwise, and T = (T,, TV), is the wind stress distribution. The subscript k + l/2 indicates variables at the interface between layers k and k + 1: wk+ l/2 is the vertical interfacial velocity; h k+1/2 = fO(+k+ 1 - &j/g;+ 1,2 represents the interface height perturbation; $k + 1,2 = ( Hk + 1$k + Hkq!tk+ l)/(Hk + Hk+ ,) is the streamfunction at the interface and g;,,,, = g (&+i - pk)/& is the reduced gravity at the interface. Bottom topography forces nonzero vertical velocity at the bottom boundary according to wN=
l/2
=
-J(h,,
!bN).
(3)
The vertical velocity at the top boundary (w,,,) is zero, consistent with the rigid lid approximation, the surface wind stress curl forcing being included as a body force over the surface layer. The vertical placement of variables is shown in Fig. 1, the model parameters in Table 1. 2.2. The fraternal-twin experiments The quasi-geostrophic model of a 3840 km square ocean basin was spun up from rest using an antisymmetric wind stress curl distribution (T,, =O; T~=T~ cos ry/L, -L
128
M.M. Rienecker, D. Adamec /Journal of Marine Systems 6 (1995) 125-143
lated altimeter data. The wind forcing, stratification and frictional parameters for the two oceans differ by 10%; other parameters such as layer depth, model resolution (Table 1) and bathymetry (Fig. 2) remain unchanged. The 10% factor was chosen merely to make the model and reference oceans measurably different; but it is also a lower limit on observational errors. The surface layer of the reference ocean was sampled as by an altimeter with orbit inclination of 45”, repeat period of 17 days and track spacing of about 130 km. The simulated data acquired in this manner, with along-track spacing of about 21.2 km, was assimilated into the QG model once per day, with all observations assumed simultaneous. The initial conditions for the assimilation cycle, conducted for 150 days, were obtained from the independent model ocean simulation. An (albeit small) ensemble of experiments was run by choosing statistically independent initial conditions from the model ocean simulation and statistically independent data from the reference ocean simulation, i.e., the initial conditions were chosen from the 3-year model integration at 30day intervals and similarly for the simulated data gathered from the reference ocean, different data being used for each assimilation experiment. A total of 32 independent experiments were conducted for the ensemble. The mean streamfunction and its root mean square (rms) variability for each layer of the model and reference oceans (without any assimilation)
4H.
_______ _____ [_.__‘________
WI/Z ___________.____________
A
W
---_____-________ ----___________-----____________ t *
W
Deformation radii Lateral Friction, A, Bottom Friction, E Wind Stress, 70 Grid spacing Time Step
Model ocean Reference ocean Model ocean Reference ocean
Model Ocean Reference Ocean
2 l/2
H3 ----________.__---____.________ v
W
3 l/2
Fig. 1. The vertical stencil of the 3-layer model.
are shown in Fig. 2. The mean streamfunction shows evidence of mesoscale variability because the statistics are based only on three years’ integration with the data saved daily. The rms variability shows the influence of the bottom topography with the local maximum in all layers (but especially layer 3) nearly 1800 km offshore in mid-basin, in addition to the expected high variability along the mid-latitude jet and along the western boundary where the jet meanders about its separation point. The amplitude of the mean and rms variability as a function of dynamical
Table 1 Model parameters Layer depths Stratification
lvz
H,=300m;H2=900m; H,=38OOm g;,, = 0.015; g;,* = 0.005 &/a = 0.017; gj,, = 0.006 R, = 30.4 km; R, = 18.7 km R, = 32.9 km; R, = 20.2 km -8x10gm4s-’ 1 x lo-‘s-1 1.0 x 10-l N m-’ 0.9 x 10-l N m-* 15 km 2 hours
M.M. Rienecker, D. Adamec /Journal
mode for each integration is shown in Fig. 3. The mid-latitude jet is strongly baroclinic with the barotropic mode confined to smaller scales in the mean. The variability influenced by the bottom topography is most apparent in the barotropic mode, consistent with the diminuition of the baroclinic modal profiles with depth. The greater thickness of layer 3 compared with the other layers leads to a larger influence of the bottom layer on the estimate of the barotropic mode, as evidenced by a comparison of Figs. 2 and 3. For the model ocean, about 30.8% of the total variance is contained in the barotropic mode, 61.2% in the first baroclinic mode and 8.0% in the second baroclinic mode. For the reference ocean, the distributions are 30.4%, 60.4% and 9.2%, respectively. The empirical orthogonal functions of streamfunction variability which are used to extend the
of Marine Systems 6 (1995) 125-143
129
surface height information through the water column are calculated from the time series of streamfunction anomaly from the mean at each grid point, i.e., the eofs are calculated as eigenvectors of the covariance matrix with (i, j)-th element < $~I,!J,! > here the < > indicates an ensemble average over time and space and the prime denotes anomaly from the temporal mean. The eofs were calculated independent of position, though to include spatially varying eofs would be a trivial enhancement of the procedure. To ensure that statistically independent data were used for the calculation, the data were subsampled in both space and time, to a 60 km horizontal grid on a monthly basis. The first eof (Fig. 4) accounts for almost 93% of the total variance throughout the water column and explains most of the variance in the upper two layers. The second eof explains just 7% of the total variance,
layer 1
layer 2
layer 3
Fig. 2. Statistics of the model ocean and the reference ocean for each layer. For the mean streamfunction contours, the dashed line indicates negative values; the contour interval is 10“ m*ss ’ m la yer 1 and 5 X lo3 rn’sC’ in layers 2 and 3. The contour interval for the streamfunction standard deviation is lo4 rn*s-l m layer 1, 5 X lo3 m2sC1 in layer 2, and 2 X lo3 m2s-’ in layer 3. The bathymetry is shown with a contour interval of 500 m. The thickened solid line is the 5000 m contour; depths shallower than 5000 m are dashed.
M.M. Rienecker, D. Adamec /Journai
130
but explains most of the variance in the bottom layer (see Table 2).
of Marine Systems 6 (1995) 125-143 Table 2 Percent variance position
W(n)
= @j(n) +ef(n),
where the correction ef(n> is derived from the observations at t,, @%z), and the model forecast field: ef(n)
=K,[W(n)
-Ha@(n)],
by EOFs
and their modal
decom-
EOF
Layer 1
Layer 2
Layer 3
Total
Mode 0
Mode 1
Mode 2
1 2 3
99.0 1.0 0.0
84.2 14.8 1.0
36.8 62.0 1.2
92.7 7.1 0.2
28.2 65.8 6.0
64.8 15.4 19.8
7.0 18.8 74.2
2.3. The assimilation scheme In general, an 01 analysis field, W(n) (denoting a vector at time t, with components c&(t,) where the subscripts i, j, k indicate the gridpoint (xi, yj, z,)), is the linear, least-squares best estimate of a prognostic field formed from a correction to the dynamical model forecast field @(n>:
explained
where K, =PH;(H,PH;+R)-I. H, defines the observational scheme at time t,; the superscript T denotes the matrix transpose; P is the forecast error covariance matrix and R is the observational error covariance. The error covariance of the analysis field, P,“, is P,” = (I-K,H,)P,
mode 0
Fig. 3. Statistics of the model ocean and the reference ocean for each mode. For the mean modal amplitude contours, the dashed line indicates negative values and the contour interval is 5 x 10’ m ’ s -’ for all modes. The contour interval for the modal amplitude standard deviation is 2.5 X lo3 m’s_’ for modes 0 and 1, and lo3 m2sC1 for mode 2.
M.M. Rienecker, D. Adamec /Journal
0.’
”
v
’
1000 -
x JI Ii $
2000.’ 3000-
*’
,’
.=’
3 4000-
5000 Fig. 4. The vertical eofs of streamfunction mined from the model ocean.
variability
deter-
where I is the identity matrix. The forecast error and the observational error are assumed to be uncorrelated. Details of this procedure are given by Miller (1986). In this study, the 01 was performed on a 2-dimensional field, so that au(n) is comprised of +t(t,) (elements on a 2-dimensional horizontal grid), and the covariance matrices are 2-dimensional. Two different methods to assimilate the surface altimeter height are considered here. For each, the upper layer streamfunction is used as a proxy for the surface altimeter height, assuming the relation: h altimeter
=
**f/g.
Throughout, the same assumptions were made with respect to the error covariances, viz., the measurement error was assumed to be white there is no correlation between measurement errors at different locations and no correlation between measurement errors and forecast errors. In addition, the model forecast error covariance, P was chosen to have an isotropic, homogeneous functional form similar to the observed streamfunction covariance:
P(r) = d[ 1 - (+)“I with a = 140 km and
e-(r’@, b =
(4)
120 km. The decorrela-
of Marine Systems 6 (1995) 125-143
131
tion scales used for this covariance are commensurate with those of the model streamfunction anomaly and assumes that the dominant contribution to the error in the forecast streamfunction comes from inadequate knowledge of the intensity and phase of the mesoscale variability. Phillips (1986) points out that this would be the case for an excellent forecast model that assimilates very little data. Certainly the cross-track updates fit into this scenario since (except in the region of the mid-latitude jet) the baroclinic information propagates across only a few grid points between updates and the update interval for a single track is about the same as the decorrelation time for the mesoscale eddy field (20-30 days). The assumed error covariance correlation scales are also commensurate with the cross-track spacing of the satellite orbits, rather than the along-track observation spacing, and is consistent with the findings of White et al. (1990b) and Rienecker and Miller (1990, i.e., if an isotropic error covariance function is used, the relevant decorrelation scale is no smaller than the larger of the meridional or zonal observation spacing. The reason for this, as Phillips (1986) points out, is that the error covariante scales estimated from observations are limited by the distance between the observations there is no information on finer scales. Only observations within a radius of 120 km were allowed to influence the analysis at any grid point since streamfunction data separated by 120 km have a correlation of only 0.1. This localization, the use of a single type of observation and the assumption of a spatially homogeneous observation noise level at a specified fraction of the model error variance mean that u2 can be an arbitrary normalization factor. It would be trivial to include a temporal variation for u2 which includes an assumed model error growth rate and the effect of reduced model forecast error variance at the points which have been corrected recently through assimilation (e.g., Rienecker and Miller, 1991). However, this has not been done here. The observational strategy is such that H,, is a matrix with non-zero (unit) entries only at model grid points which lie beneath a satellite orbit (see White et al., 1990a for a typical sampling strategy on the model grid).
132
M.M. Rienecker, D. Adamec /Journal
It should be pointed out that the use of approximate error covariance functions reduces the method to being a statistical interpolation rather than a statistically optimal interpolation. Daley (1991) discusses the impact of errors in assumed covariance functions. In practice, the use of assumed funtional forms for error covariances is usually the best that can be done and is still found to be effective. The two different techniques considered here assimilate: (1) the estimated amplitude, (Y~(x~,Yj, t,), of the dominant vertical eof of streamfunction variability, so that ~ij(t,> = (~i(x~, Yj, t,); and 2) the surface streamfunction, so that ~$~~(t,)= @i(xj, Yjy t,). Th e f irs t is referred to hereafter as the eof-based method, the second as the surfacebased method.
of Marine Systems 6 (1995) 125-143
The eof-based method
The assimilation cycle estimates (bG(t,) = c#xj, .Yj, f,J for each observation by assuming that all of the observation signal is contained in the first eof: (y;)(xi? Yj? ‘n) = Ilrf(‘i? Yj? ‘n)/511,
(5)
where tlk is the k-th component of the I-th eof of streamfunction variability and (cl$xi, yj, t,) is the surface layer streamfunction (surface height) observation sampled from the reference ocean. Since the eofs form a complete set of basis functions for the 3-layer streamfunction, the amplitude of the first eof from the model [c##,)] is obtained by projection: iv=3 af(XiY
Yj7 ‘n) =
C
(clk(xi,
k=l
Yj, fn)51k.
(6)
DAY 5000 Reference
Difference
layer 1
layer
2
layer 3
Fig. 5. The initial assimilation streamfunction analyses, the streamfunction from the reference ocean and the difference between the two, for each layer from one of the experiments conducted. The contour interval is 2 X lot, lo4 and lo4 m*s-’ in layers 1, 2 and 3, respectively. Contours less (greater) than - (+) 4 X lo4 in layer 1, and - (+) 2 X lo4 in layers 2 and 3 are shaded blue (orange).
M.M. Rienecker, D. Adamec /Journal
two-gyre experiment, but also marked differences in the meandering pattern of the mid-latitude jet and the intensity and placement of the eddies. The differences extend throughout the water column with similar scales (up to about 700 km). The inital root mean square difference (error> between the analysis and the reference field is 64% of the streamfunction variance of the reference ocean in the upper layer, 108% in the second layer and 89% in the bottom layer (the nrms statistic in Table 3). The seemingly high correlation (0.81 in the surface layer, 0.45 in the second layer and 0.58 in the bottom layer) show the similarities expected in the streamfunction distribution associated with the time-invariant surface forcing specified (only the basin-wide mean has been removed), as well as the influence of the bottom topography in steering the flow. After a complete orbit cycle (day 5020, Figs. 6a, 7a), both the eof method and the surface method have produced analyses much closer to the reference ocean, with differences smaller in scale (now about 200 km) throughout the water column for the eof method. The differences have been noticeably reduced in the far field of the mid-latitude jet (as also found by White et al., 1990a). The global rms differences have decreased by 52%, 41% and 11% of their initial values in each of the three layers, respectively, for the eof-based analyses (Table 3). In contrast, the rms error has decreased by 38% for the surface layer, but increased slightly in the subsurface layers for the surface-based analyses, and the differences are slightly larger scale. After three
After the 01 analysis is produced, the model integration proceeds with I,!J~ where
N=3 +
C af(Xi, 1=2
Yj9 fn)5[k,
133
of Marine Systems 6 (1995) 125-143
(7)
with (Y{,1= 2,3, determined analogous to Eq. (6). A variation on this would partition the observation signal across all eofs according to the variance explained by each, as in Dombrowsky and De Mey (1992); however since the higher eofs explain less than the assumed noise level (10% of the signal variance), such an approach is not warranted here. The inclusion of higher eofs would undoubtedly require a different form, or different parameters, for P in Eq. (4). The surface-based method
For this method, the 01 operates directly on the surface streamfuntion, $i(x,, yj, t,>. After the 01, the model integration proceeds with the analysis streamfunction in the surface layer, but the subsurface streamfunction unchanged.
3. Results
A comparison of the initial field used for one of the assimilation experiments (beginning model day 5000, Fig. 5) with the corresponding reference field from the series used for simulated data shows the similarities expected of this type of
Table 3 Root mean square error and correlation measures for one experiment. The root mean square difference between the analyses and the reference ocean is normalized by the global variance in the reference ocean (nrms) and given as a percentage Day
0 20 50 100
nrms C nrms C nrms C nrms C
EOF method
Surface method
Layer 1
Layer 2
Layer 3
Layer 1
Layer 2
Layer 3
64 0.81 31 0.96 32 0.96 28 0.97
108 0.45 64 0.84 71 0.81 51 0.87
89 0.58 75 0.75 80 0.76 55 0.86
64 0.81 40 0.92 44 0.90 58 0.83
108 0.45 115 0.48 141 0.54 239 0.26
89 0.58 101 0.53 132 0.62 250 0.45
MM. Rienecker, D. Adamec /Journal
134
of Marine Systems 6 (1995) 125-143
DAY 5020 Difference
layer
1
layer
2
layer 3
DAY 5050 Difference
layer 1
layer 2
layer 3
(a)
M.M. Rienecker, D. Adamec /Journal
of Marine Systems 6 (1995) 125-143
135
Difference
layer 1
layer 2
layer 3
Fig. 6. The analyses from the eof-based assimilation technique, streamfunction from the reference ocean and the difference between the two, for each layer. The fields after (a) 20 days (b) 50 days and (c) 100 days of assimilation. Contour intervals and shading as for Fig. 5.
repeat cycles (day 5050), it appears that the similarity between the eof-based analysis and the reference ocean has not markedly improved, except that the region of greatest differences is more tightly bound to the jet (Fig. 6b). In contrast, the differences in the surface-based analysis have increased in scale to about 1000 km in the subsurface layers and the magnitude of the differences has increased substantially, to (globally) 141% of the natural variability for layer 2 and 132% for layer 3 (Fig. 7b, Table 3). After 100 days of assimilation, the eof-based analyses are visually in very close agreement with the reference ocean in all layers (Fig. 6c) and the global statistics confirm the good agreement in all layers (Table 3), with a reduction of the initial error by 56% in the surface layer, 53% in the second layer
and 38% in the bottom layer. Whereas there are some features in the surface-based analysis which compare well with the reference field, there are also large discrepancies (Fig. 7c, note the difference in contour interval) and the correlation between the analyses and the reference field has decreased significantly (Table 3). In the surface layer, the intense discrepancies are still small scale, but in the subsurface layers the differences are larger scale (up to 1000 km) and are not confined to the region of the jet. This is evident in the comparisons of Fig. 8 where the absolute error of the analysis field is normalized by the standard deviation of the variability at each grid point in the reference ocean. After 100 days of assimilation, the improvement in the eof-based analysis is global, whereas the disagreement (note
M.M. Rienecker, D. Adamec/Journal
136
of Marine Systems 6 (1995) 125-143
DAY 5020
layer
1
layer
2
layer
3
DAY 5050 Reference
layer 1
layer
2
layer
3
Difference
M. M. Rienecker, D. Adamec /Journal
of Marine Systems 6 (1995) 125-143
137
DAY 5100 Difference
layer 1
layer 2
layer 3
Fig. 7. The analyses from the surface-based assimilation technique, streamfunction from the reference ocean and the difference between the two, for each layer. The fields after (a) 20 days (b) 50 days and (c) 100 days of assimilation. The contour intervals are as for Fig. 5, except for 100 days when the contour intetval in layers 2 and 3 are the same as for layer 1.
again the difference in contour interval for the surface- based method at day 5100) in the subsurface layers of the surface-based method extend throughout a large fraction of the domain, but are most intense along the western boundary. The large relative amplitude along the eastern boundary is merely due to the low variance (normalizing factor) there, as seen in Fig. 3. In this region. linear dynamics apply and the application of a correction to the forecast streamfunction is akin to wind stress curl forcing of Rossby waves which are then clearly seen in the surface signal. The time series of the root mean square difference between the analyses and the reference ocean for these experiments are shown in Fig. 9a. Most of the reduction in error occurs very quickly, over the first satellite repeat cycle. After two
repeat cycles (34 days), the error level for the eof-based analyses is fairly constant, with excursions about the asymptote due to mesoscale activity not resolved by the satellite observations. Finer observation spacing would lower the error bound, as shown by Holland and Malanotte-Rizzoli (1989). The error level in the bottom layer slowly decreases in time, and it would likely take much longer (gyre recirculation timescales) for the error level to reach its minimum. In contrast, the error levels for the surface-based method decrease dramatically only in the surface layer. The errors stay low for about 50 days in all layers, but then start to increase without any apparent bound. Eventually, as found by White et al (1990a), the initial model field would have given a closer agreement with the reference ocean than the
MM
138
Rienecker, D. Adamec /Journal
qf Marine Systems 6 (1995) 125-143
Day 5020
EOF Method
Surface
Method
EOF Method
Day 5100 Surface
Method
layer 1
layer
2
layer
3
Fig. 8. The absolute difference between the analyses and the reference ocean, normalized by the standard deviation of the streamfunction variability at each point in the reference ocean. The contour interval is 2, except for the last column where the contour interval is 4. The shading is blue for values between 2 and 4 (4 and 8 for the last column) and orange for values greater than 4 (greater than 8 for the last column).
analysis field and the model should be re-initialized. The results for this single experiment is typical of that for each experiment in the ensemble (Fig. 9b), indicating the results are robust, at least for the dynamical scenario (wind forcing, bottom topography, stratification, resolution) considered here. The rrns error statistics for the dynamical modes (Fig. 912,d) explain the differences in behaviour of the two methods. The eof-based method shows a marked decrease in the error in the barotropic and first baroclinic modes and a slight decrease for the second baroclinic mode
over the entire assimilation experiment. In contrast, the surface-based method shows initial improvement in the barotropic and first baroclinic modes, but subsequently, the error is reduced only in the first baroclinic mode for which the error remains low for about 100 days. The error in the barotropic mode increases markedly after one or two satellite repeat cycles and the error in the second baroclinic mode increases throughout the entire assimilation experiment. The explanation for the differences in modal errors lies in the decomposition of the correction, ef(n), to the model forecast streamfunction. For
M.M. Rienecker, D. Adamec /Journal
of Marine Systems 6 (1995) 125-143
139 (b)
.-’ .A____ r”....“‘.-...........~..................,...................~. _____________________
5.0-103-~L.,..
50
100 assimilation day
0: 0
150
I 50
I 100 assimilation day
7
150
k.4
ot. 0
1.
50
I”‘.1
,‘,
150
100 assimilation day
ot 0
.
I
.
.
.
.
50
.
100
.
...1 150
assimilation day
(4
:..................................................,,..................... ===z
~____-_---r=========;; !
0
50
,...
I
100 assimilation day
I.
150
Fig. 9. Time series of root mean square error during the assimilation cycle. (a) A single assimilation experiment, the solid line is the streamfunction error in layer 1, the dashed line is the error in layer 2 and the dotted line is the error in layer 3. The thicker lines are for the eof-based method, the thinner lines for the surface method. (b) As for (a), but the ensemble average error. (c) As for (a), but for the modes. The solid line is for the barotropic mode, the dashed line for the first baroclinic mode, and the dashed line for the second baroclinic mode. Cd) As for Cc) but for the ensemble average error. (e) As for (a) but for potential vorticity (scaled by 1O6).
M.M. Rienecker, D. Adamec /Journal
140
the eof-based method, 28.2% of the variance goes to the barotropic mode, 64.8% to the first baroclinic mode and 7.0% to the second baroclinic mode, very close to the distribution for the reference ocean. In contrast, for the surface-based method, only 6% of the variance goes to the barotropic mode, 65% to the first baroclinic mode and 29% to the second baroclinic mode. The continuous addition of anomalous energy in the higher baroclinic mode apparently leads to larger transfers to barotropic energy through the various energy cascades, some of which may be enhanced by topography (e.g. Adamec, 1988). The larger subsurface scales in the difference fields for the surface-based method (Fig. 7) are associated with the error in the barotropic mode and its scales as set by topography. For the surface layer, the larger scales are well-defined by the observation strategy (viz., scales of about 300 km for the Geosat track spacing, as shown by White et al., 1990b), so the errors remain in the smaller (unresolved) scales. It is of interest to compare the present vertical projection technique with that used by Haines (1991) wherein the streamfunction in the surface layer is replaced with the observed streamfunction and the subsurface streamfunction is adjusted so that the subsurface potential vorticity, &, is unchanged, where
fll
!L= v*V+,+f(y)
- #Q-i,*
-&Cc,,*)
k +s
f&B i’.‘,k HN .
This process entails solving elliptic equations for each subsurface streamfunction correction: V26$, - (Y& + Y22,,)& + Y&V3 = -Y;,,~&~ V2S6, -
Y32,2W3
+
d,*w2
=
0,
where, using Haines’ notation, yi k +, = ft/g; +1,2Hk and c?$~ is the streamfunction correction in layer k. With a spectral (horizontal) decomposition and the model parameters used in this study, it is easy to show that for very small wavenumbers, the response is S& = ?j+2 = S+i, i.e., the constant potential vorticity method projects the surface signal onto the
of Marine Systems 6 (1995) 125-143
barotropic mode. For wavenumbers = (70 km)-‘, the response is 8$s = 0.2564, and SI&*2: 0.46$,, ratios which are very similar to those for the first eof in this study which has 5,s = 0.205,, and 5i2 = 0.375ii. For smaller scales, the surface signal has much reduced influence on the subsurface layers and at the first baroclinic radius of deformation, the influence is very small, particularly in the bottom layer: S$, = O.O6S$, and &,!I, = 0.2s*,. In an identical-twin, flat bottom ocean study, Haines found that the error in the bottom streamfunction increased during assimilation but that afterwards the bottom Ekman drag and the interfacial eddy stresses force the bottom layer towards the control ocean. The damping effect of bottom drag and the “correcting influence of the currents generated from above” help redistribute the potential vorticity in the bottom layer. However, the dynamical coupling is two-way; i.e., the variability (errors) in the bottom layer also force variability in the middle layer(s) (see Lozier and Riser, 1990, for the long-term balances in eddy driving of subsurface layers in a flat-bottom ocean) and the errors may magnify in time. In addition, because of the thickness of the bottom layer, errors in potential vorticity in this layer, either in the initial conditions or introduced through assimilation, have predominant influence on the barotropic mode. This mode has the strongest interaction with the bottom topography. Outside of smooth, large-scale bottom slopes, topographic interactions tend to change the energy cascade toward smaller scale barotropic flow as well as surface-intensified baroclinic flow (e.g. Adamec, 1988), the latter requiring further adjustments in the upper layer flow in the assimilation process. For the eof method, the surface information is projected onto the subsurface layers according to the long-term statistical relationship between the layers - these statistics automatically encompass the energy cascades consistent with the underlying dynamics, including the surface and bottom interactions, and any statistically significant potential vorticity homogenization. Of course, the streamfunction ratios are not correct at all places and time. For both the constant potential vortic-
M.M. Rienecker, D. Adamec /Journal
ity method and the eof method, the reliance is on the model to adjust the incorrect potential vorticity according to the new analysis streamfunction in each layer and the dynamical balances of the model. The difference between the potential vorticity for the analyses and the reference ocean is shown in Fig. 9e. For this case study, the potential vorticity difference in the lower layers remains fairly constant for the eof-based method. The assimilation procedure reduces the potential vorticity error in the upper layer by 39% after 150 days and increases the potential vorticity error by 15% in the middle layer and by 0.4% in the bottom layer. Maps of potential vorticity (not shown) show that the large pool of fairly homogeneous potential vorticity in the middle layer, beneath the surface mid-latitude jet, is not destroyed by the eof-based analyses. The bottom layer has no such homogeneous region because of the bottom topography.
4. Discussion It has been demonstrated that optimal interpolation and the use of a single model-based vertical eof of streamfunction variability are effective and efficient for assimilating surface altimeter data into fraternal-twin simulations of a 3-layer quasigeostrophic model with topography. After 100 days of assimilation (about 6 satellite repeat cycles), the error in the upper two layers has been reduced by over 50% and that in the bottom layer by 38%. The essence of the method is that the eofs capture the statistics of the dynamical balances in the model and ensure that this balance is not inappropriately disturbed during the assimilation process. The analysis fields have been filtered according to the model’s dynamical and statistical balances. In contrast, a method which assimilates the information at the surface only and relies on the model to extend the information throughout the water column imposes an incorrect distribution of energy on the fundamental dynamical modes of the model; the reliance then is on the model dynamics to redistribute the energy according to its own nonlinear cascade. Of course, if the assimilation were curtailed, the
of Marine Systems 6 (1995) 125-143
141
model would do exactly that, reproducing, after a few eddy turn-over times, an energy distribution more like that in the model ocean. This suggests that there may be some benefit (A. Robinson, pers. commun.) in assimilating data infrequently and allowing the model to redistribute the energy according to its own dynamical balances. For the eof-based method, such considerations are not necessary and the data may be used close to its observation time to try to preserve the phase information as closely as possible. The eof-based method effectively reduces the degrees of freedom and so the computational burden of the assimilation process. The success of using a single eof for the entire domain relies somewhat on the fact that here this single eof explains more than 90% of the entire lateral and vertical streamfunction variance. Such high levels of variance explained by a single mode are not unprecedented in data: observations in limited area domains (e.g., Rienecker et al., 1987; Dombrowsky and De Mey, 1992) show dominance by a single local vertical mode of dynamic topography. In practice, the implementation of additional and spatially varying eofs would pose no difficulties to the use of this technique for effective assimilation of surface information. The implementation requires the specification of a horizontal forecast error covariance function for each eof. The underlying assumption is that the 3-dimensional error covariante is separable (horizontal from vertical) so that the analysis may be decoupled in a similar manner to the vertical decoupling of the dynamical equations by projection onto dynamical vertical modes. The use of eofs (statistical modes) is more useful than dynamical modes if (as in this study) the variance is predominantly in a few statistical modes. For special forms of the vertical observation error covariance matrix (see Daley, 1991, Ch. 41, it can be shown that the eofs of the background error covariance (here assumed to be the same as the eofs of streamfunction variability) are also eofs (and so an efficient characterization) of the analysis error covariance. The full Kalman filter would evolve this analysis error covariance according to the nonlinear, vertically-coupled dynamics. However, as long as the system (model) error is not dominant, it might be expected that
142
M.M. Rienecker, D. Adamec /Journal
(unless a lot of data is assimilated often) the initial eofs would remain a good representation of the error covariance as they are of the streamfunction variability which has evolved under the same dynamics. Even with a model with much higher vertical resolution, it is unlikely that a large number of eofs would be required. The use of a small number of eofs of the forecast error covariance matrix has been found efficacious in the National Meteorological Center’s global data assimilation system (Parrish and Derber, 1992). For the ocean, in a joint eof decomposition of fields of temperature, salinity, oxygen and nutrient observations in the North Atlantic, Fukumori and Wunsch (1991) were able to explain 92% of the variance with 6 modes, 86% with 4 modes and 68% with only 2 modes. For model-based eofs, the explained variance for the low modes is likely to be higher than for observations because of lower noise levels. Lack of knowledge of surface forcing parameters and uncertainty in fundamental model parameters could lead to a model-estimated eofs which are substantially different from that supported by analysis fields influenced by a lot of data (because of different ratios of barotropic to baroclinic energy etc). For this case, new eofs could be reestimated after many assimilation cycles. In the assimilation of satellite observations (rather than simulated data) it will be important to include along-track observational error covariances.
Acknowledgements The authors gratefully acknowledge enlightening discussions with and suggestions from Mike Bell, Bob Haney, Robert Miller, Allan Robinson, and John Marshall. The comments of two anonymous reviewers helped clarify the presentation and are greatly appreciated. This research was supported by the Ocean Modeling and Data Analysis Program under NASA’s Climate and Hydrologic Systems Branch. Computing was carried out at the NASA Center for Computational Sciences at Goddard Space Flight Center.
of Marine Systems 6 (1995) 125-143
References Adamec, D., 1988. Numerical simulations of the effects of seamounts and vertical resolution on strong ocean flows. J. Phys. Oceanogr., 18: 258-269. Bengtsson, L., Ghil, M. and Kallen, E. (Editors), 1981. Dynamic Meteorology: Data Assimilation Methods. Springer, New York, 330 pp. Cummins, P. and Mysak, L.A., 1988. A Quasi-geostrophic circulation model of the Northeast Pacific, Part I: A preliminary numerical experiment, J. Phys. Oceanogr., 18: 1261-1286. Daley, R., 1991. Atmospheric Data Analysis. Cambridge University Press, 457 pp. De Mey, P. and Robinson, A.R., 1987. Assimilation of altimeter eddy fields in a limited-area quasigeostrophic model. J. Phys. Oceanogr., 17: 2280-2293. Dombrowsky, E. and De Mey, P., 1992. Continuous assimilation in an open domain of the Northeast Atlantic 1. Methodology and application to AthenaA-88. J. Geophys. Res., 97: 9719-9731. Fukumori, I. and Wunsch, C., 1991. Efficient representation of the North Atlantic hydrographic and chemical distributions. Prog. Oceanogr., 27: 111-195. Haines, K., 1991. A direct method for assimilating sea surface height data into ocean models with adjustments to the deep circulation, J. Phys. Oceanogr., 21: 843-868. Holland, W.R., 1989. Altimeter data assimilation into ocean circulation models - some preliminary results. In: D.L.T. Anderson and J. Willebrand (Editors), Oceanic Circulation Models: Combining Data and Dynamics. Kluwer, Dordrechts, 605 pp. Holland, W.R. and Malanotte-Rizzoli, P., 1989. Assimilation of altimeter data into an ocean circulation model: space versus time resolution studies. J. Phys. Oceanogr., 19: 1507-1534. Hurlburt, H., 1986. Dynamic transfer of simulated altimeter data into subsurface information by a numerical ocean model. J. Geophys. Res., 91: 2372-2400. Lozier, MS. and Riser, S., 1990. Potential vorticity sources and sinks in a quasi-geostrophic ocean: beyond Western Boundary Currents. J. Phys. Oceanogr., 20: 1608-1627. Marshall, J. and Berry, P., 1991. Ocean modelling studies inspired by the Seasat altimeter. Int. J. Remote Sensing, 12: 1631-1638. Mellor, G.L. and Ezer, T., 1991. A Gulf Stream model and an Altimetry assimilation scheme. J. Geophys. Res., 96: 8779-8795. Miller, R.N., 1986. Toward the application of the Kalman filter to regional open ocean modeling. J. Phys. Oceanogr., 16: 72-86. Parrish, D.F. and Derber, J.C, 1992. The National Meteorological Center’s Spectral Statistical-Interpolation Analysis Scheme. Mon. Weather Rev., 120: 1747-1763. Phillips, N.A., 1986. The spatial statistics of random geostrophic modes and first-guess errors. Tellus, 38A: 314-332.
M.M. Rienecker, D. Adamec /Journal
Rienecker, M.M. and Miller, R.N., 1991. Ocean data assimilation using Optimal Interpolation with a quasi-geostrophic model. J. Geophys. Res., 96, l&093-15,103. Rienecker, M.M., Mooers, C.N.K. and Robinson, A.R., 1987. Dynamical interpolation and forecast of the evolution of mesoscale features off Northern California. J. Phys. Oceanogr., 17: 1189-1213. Verron, J., 1990. Altimeter data assimilation into an ocean model: sensitivity to orbital parameters. J. Geophys. Res., 95, 11,443-11,459. Verron, J., 1992. Nudging satellite altimeter data into quasi-
of Marine Sjstems 6 (1995) 125-143
143
geostrophic ocean models. J. Geophys. Res., 97: 74797491. White, W.B., Tai, C.-K. and Holland, W.R., 1990a. Continuous assimilation of simulated Geosat altimetric sea level into an eddy-resolving numerical ocean model, 1, Sea level differences. J. Geophys. Res., 95: 3219-3234. White, W.B., Tai, C.-K. and Holland, W.R., 1990b. Ccmtinuous assimilation of simulated Geosat altimetric sea level into an eddy-resolving numerical ocean model, 2, Referenced Sea level differences. J. Geophys. Res., 95: 32353251.
ELSEVIER
OF
MARINE SYSTEMS Journal of Marine Systems 6 (1995) 125-143
Assimilation of altimeter data into a quasigeostrophic model using optimal interpolation and eofs
ocean
M.M. Rienecker, D. Adamec Code 971, Laboratory for Hydrospheric Processes, NASA /Goddard
Space Flight Center, Greenbelt, MD 20771, USA
Received 20 September, 1993; revised and accepted 24 March, 1994
Abstract An ensemble of fraternal-twin experiments is used to assess the utility of optimal interpolation and model-based vertical empirical orthogonal functions (eofs) of streamfunction variability to assimilate satellite altimeter data into ocean models. Simulated altimeter data are assimilated into a basin-wide 3-layer quasi-geostrophic model with a horizontal grid spacing of 15 km. The effects of bottom topography are included and the model is forced by a wind stress curl distribution which is constant in time. The simulated data are extracted, along altimeter tracks with spatial and temporal characteristics of Geosat, from a reference model ocean with a slightly different climatology from that generated by the model used for assimilation. The use of vertical eofs determined from the model-generated streamfunction variability is shown to be effective in aiding the model’s dynamical extrapolation of the surface information throughout the rest of the water column. After a single repeat cycle (17 days), the analysis errors are reduced markedly from the initial level, by 52% in the surface layer, 41% in the second layer and 11% in the bottom layer. The largest differences between the assimilation analysis and the reference ocean are found in the nonlinear regime of the mid-latitude jet in all layers. After 100 days of assimilation, the error in the upper two layers has been reduced by over 50% and that in the bottom layer by 38%. The essence of the method is that the eofs capture the statistics of the dynamical balances in the model and ensure that this balance is not inappropriately disturbed during the assimilation process. This statistical balance includes any potential vorticity homogeneity which may be associated with the eddy stirring by mid-latitude surface jets.
1. Introduction The assimilation of oceanographic data into numerical models is essential not only for predictive purposes, but also for diagnostic purposes, to help further the understanding of the ocean circulation on regional, basin and global scales. Given the sparsity and asynopticity of hydrographic and current meter data, oceanographers look forward to the availability of more numerous data through remote sensing. Satellites provide many measurements of several ocean variables;
but the measurement of the ocean surface topography by satellite altimeters provides the most direct connection to ocean dynamics since the altimeter signal is essentially a measure of vertically integrated mass. The problem arises as to how best to use this surface measurement to gain information about subsurface ocean variability. This paper describes a technique for projecting the surface information provided by the satellite altimeter to the ocean subsurface in an efficient manner and compares the results with a similar technique which relies entirely on the dynamical
0924-7963/9.5/$09.50 0 1995 Elsevier Science B.V. All rights reserved SSDI 0924-7963(94)00013-2
126
MM. Rienecker, D. Adamec /Journal
model to derive the subsurface variability consistent with that at the surface. Altimeter data has been incorporated into dynamical models by either assimilating the information only at the surface of the model (e.g., Hurlburt, 1986; Holland, 1989; Holland and Malanotte-Rizzoli, 1989; Verron, 1990, 1992; White et al, 1990a, b; Marshall and Berry, 1991) or by projecting the information in some manner in the vertical (e.g., De Mey and Robinson, 1987; Haines, 1991; Mellor and Ezer, 1991). The projection techniques have used either statistical information or dynamical “constraints”. The former methods include the use of empirical orthogonal functions of density perturbations determined from observations (De Mey and Robinson, 1987) or the use of regression techniques to estimate subsurface patterns (Mellor and Ezer, 1991). Holland and Malanotte-Rizzoli (1989) rely on model dynamics to project the surface information throughout the water column, viz., they assimilate along-track potential vorticity (i.e., the contribution to the potential vorticity from along-track derivatives of surface height) in the surface layer of their quasi-geostrophic (QG) model. Since the QG layer streamfunction is coupled vertically through the definition of QG potential vorticity, the surface information is passed immediately through the water column. In contrast, Haines (1991) inserts the surface streamfunction in the topmost layer only of the QG model and applies constraints on the potential vorticity in the subsurface layers and finds that the surface information is effectively projected through most of the layers. Different techniques are also used to incorporate the horizontal distribution of information. The most common methods are direct insertion, nudging and optimal interpolation. For direct insertion, the model streamfunction is simply replaced by the observation, as, e.g. in Haines (1991) who assumes that observations are available at all model grid points. For nudging an extra forcing term, proportional to the difference between the model variable and the observation, is added to the dynamical equations (e.g. Holland and Malanotte-Rizzoli, 1989). Optimal Interpolation (011, the data assimilation technique most com-
of Marine Systems 6 (1995) 125-143
monly used in operational numerical weather prediction (see, e.g. Bengtsson et al., 19811, produces an analysis field which is a linear combination of model-generated streamfunction and observed streamfunction. By using weights based on estimated error covariances for the model forecasts and the measurements the linear combination minimizes, in a least-squares sense, the errors in the analysis. White et al. (1990a, b), in identical-twin experiments, assimilated simulated Geosat altimeter observations into a 3-layer eddy-resolving, flat bottom, square basin QG model using 01 in the surface layer of the model only. Their first study assimilated altimeter differences into a model with initial fields 1% in error from the reference fields and found that the nonlinear regime of the mid-latitude jet limited the ability of the observations to constrain the model, so that the model had to be re-initialized periodically. Their second study placed stronger controls on the influence of the observations through the error covariance and in repeated use of the observations. The simulated altimeter observations were referenced to the model mean streamfunction and the study found that the assimilation procedure was effective in increasing the correlation between analyses and the reference ocean for all layers, however, the integrations were conducted for only 50 days. Marshall and Berry (1991) note that whereas surface information appears to be useful in constraining the subsurface flow in a two-layer model (depending on the satellite repeat time), this ability is severely limited in models with higher vertical resolution. It is apparent that techniques should be sought to the aid the model in its dynamical extrapolation. In this paper, an ensemble of fraternal-twin experiments is used to assess the assimilation, by optimal interpolation, of simulated surface altimeter data into a 34ayer square-basin QG model with bottom topography, forced by a wind stress curl distribution which is constant in time. The simulated data, along altimeter tracks with spatial and temporal characteristics of Geosat, are extracted from a time-series of surface-layer streamfunction in a reference ocean generated by
M.M. Rienecker, D. Adamec /Journal
the QG model with slightly different model parameters. The results from assimilation of the projection of the data onto empirical orthogonal functions (eofs) of the model streamfunction is compared with those from assimilation of the data into the surface layer only of the model. The use of model-generated eofs can be considered a compromise between the use of eofs derived by observations as in De Mey and Robinson (1987) and Dombrowski and De Mey (1992) and the use of a specific dynamical constraint, as by Haines (1991). The resulting analysis fields are therefore filtered according to the model’s dynamics and statistics. This paper is organized as follows: the assimilation system and the fraternal-twin experiment design is described in Section 2, the results are presented in Section 3 and conclusions summarized in Section 4.
2. The assimilation
system
2.1. The QG model The dynamical model used in the assimilation system is a three-layer QG model with solution technique as described by Cummins and Mysak (1988). For an N-layer model, the QG streamfunction is found as the solution to
+
fll j-+-l,2
-
wk+1/2
)
-44%
k
curl7 - $@v2*,
+ a,,--
7
PI+,
k=l,...,N;
(I)
and ;cs,
-
&+I)
=J(vh
-
h+*7
4&+1,2)
&+1/2 -
k=l,...,N-1.
---wk+l/2,
fo
(2)
of Marine Systems 6 (1995) 125-143
121
In the above, V2 is the horizontal Laplacian operator, J is the Jacobian operator, f = f0 + fly is the Coriolis parameter on the P-plane centered at 35”N, p is the density; g is the acceleration due to gravity; h, is the bottom topography; Hk are the average layer thicknesses, I?,, = 1, if 1= k, 0 otherwise, and T = (T,, TV), is the wind stress distribution. The subscript k + l/2 indicates variables at the interface between layers k and k + 1: wk+ l/2 is the vertical interfacial velocity; h k+1/2 = fO(+k+ 1 - &j/g;+ 1,2 represents the interface height perturbation; $k + 1,2 = ( Hk + 1$k + Hkq!tk+ l)/(Hk + Hk+ ,) is the streamfunction at the interface and g;,,,, = g (&+i - pk)/& is the reduced gravity at the interface. Bottom topography forces nonzero vertical velocity at the bottom boundary according to wN=
l/2
=
-J(h,,
!bN).
(3)
The vertical velocity at the top boundary (w,,,) is zero, consistent with the rigid lid approximation, the surface wind stress curl forcing being included as a body force over the surface layer. The vertical placement of variables is shown in Fig. 1, the model parameters in Table 1. 2.2. The fraternal-twin experiments The quasi-geostrophic model of a 3840 km square ocean basin was spun up from rest using an antisymmetric wind stress curl distribution (T,, =O; T~=T~ cos ry/L, -L
128
M.M. Rienecker, D. Adamec /Journal of Marine Systems 6 (1995) 125-143
lated altimeter data. The wind forcing, stratification and frictional parameters for the two oceans differ by 10%; other parameters such as layer depth, model resolution (Table 1) and bathymetry (Fig. 2) remain unchanged. The 10% factor was chosen merely to make the model and reference oceans measurably different; but it is also a lower limit on observational errors. The surface layer of the reference ocean was sampled as by an altimeter with orbit inclination of 45”, repeat period of 17 days and track spacing of about 130 km. The simulated data acquired in this manner, with along-track spacing of about 21.2 km, was assimilated into the QG model once per day, with all observations assumed simultaneous. The initial conditions for the assimilation cycle, conducted for 150 days, were obtained from the independent model ocean simulation. An (albeit small) ensemble of experiments was run by choosing statistically independent initial conditions from the model ocean simulation and statistically independent data from the reference ocean simulation, i.e., the initial conditions were chosen from the 3-year model integration at 30day intervals and similarly for the simulated data gathered from the reference ocean, different data being used for each assimilation experiment. A total of 32 independent experiments were conducted for the ensemble. The mean streamfunction and its root mean square (rms) variability for each layer of the model and reference oceans (without any assimilation)
4H.
_______ _____ [_.__‘________
WI/Z ___________.____________
A
W
---_____-________ ----___________-----____________ t *
W
Deformation radii Lateral Friction, A, Bottom Friction, E Wind Stress, 70 Grid spacing Time Step
Model ocean Reference ocean Model ocean Reference ocean
Model Ocean Reference Ocean
2 l/2
H3 ----________.__---____.________ v
W
3 l/2
Fig. 1. The vertical stencil of the 3-layer model.
are shown in Fig. 2. The mean streamfunction shows evidence of mesoscale variability because the statistics are based only on three years’ integration with the data saved daily. The rms variability shows the influence of the bottom topography with the local maximum in all layers (but especially layer 3) nearly 1800 km offshore in mid-basin, in addition to the expected high variability along the mid-latitude jet and along the western boundary where the jet meanders about its separation point. The amplitude of the mean and rms variability as a function of dynamical
Table 1 Model parameters Layer depths Stratification
lvz
H,=300m;H2=900m; H,=38OOm g;,, = 0.015; g;,* = 0.005 &/a = 0.017; gj,, = 0.006 R, = 30.4 km; R, = 18.7 km R, = 32.9 km; R, = 20.2 km -8x10gm4s-’ 1 x lo-‘s-1 1.0 x 10-l N m-’ 0.9 x 10-l N m-* 15 km 2 hours
M.M. Rienecker, D. Adamec /Journal
mode for each integration is shown in Fig. 3. The mid-latitude jet is strongly baroclinic with the barotropic mode confined to smaller scales in the mean. The variability influenced by the bottom topography is most apparent in the barotropic mode, consistent with the diminuition of the baroclinic modal profiles with depth. The greater thickness of layer 3 compared with the other layers leads to a larger influence of the bottom layer on the estimate of the barotropic mode, as evidenced by a comparison of Figs. 2 and 3. For the model ocean, about 30.8% of the total variance is contained in the barotropic mode, 61.2% in the first baroclinic mode and 8.0% in the second baroclinic mode. For the reference ocean, the distributions are 30.4%, 60.4% and 9.2%, respectively. The empirical orthogonal functions of streamfunction variability which are used to extend the
of Marine Systems 6 (1995) 125-143
129
surface height information through the water column are calculated from the time series of streamfunction anomaly from the mean at each grid point, i.e., the eofs are calculated as eigenvectors of the covariance matrix with (i, j)-th element < $~I,!J,! > here the < > indicates an ensemble average over time and space and the prime denotes anomaly from the temporal mean. The eofs were calculated independent of position, though to include spatially varying eofs would be a trivial enhancement of the procedure. To ensure that statistically independent data were used for the calculation, the data were subsampled in both space and time, to a 60 km horizontal grid on a monthly basis. The first eof (Fig. 4) accounts for almost 93% of the total variance throughout the water column and explains most of the variance in the upper two layers. The second eof explains just 7% of the total variance,
layer 1
layer 2
layer 3
Fig. 2. Statistics of the model ocean and the reference ocean for each layer. For the mean streamfunction contours, the dashed line indicates negative values; the contour interval is 10“ m*ss ’ m la yer 1 and 5 X lo3 rn’sC’ in layers 2 and 3. The contour interval for the streamfunction standard deviation is lo4 rn*s-l m layer 1, 5 X lo3 m2sC1 in layer 2, and 2 X lo3 m2s-’ in layer 3. The bathymetry is shown with a contour interval of 500 m. The thickened solid line is the 5000 m contour; depths shallower than 5000 m are dashed.
M.M. Rienecker, D. Adamec /Journai
130
but explains most of the variance in the bottom layer (see Table 2).
of Marine Systems 6 (1995) 125-143 Table 2 Percent variance position
W(n)
= @j(n) +ef(n),
where the correction ef(n> is derived from the observations at t,, @%z), and the model forecast field: ef(n)
=K,[W(n)
-Ha@(n)],
by EOFs
and their modal
decom-
EOF
Layer 1
Layer 2
Layer 3
Total
Mode 0
Mode 1
Mode 2
1 2 3
99.0 1.0 0.0
84.2 14.8 1.0
36.8 62.0 1.2
92.7 7.1 0.2
28.2 65.8 6.0
64.8 15.4 19.8
7.0 18.8 74.2
2.3. The assimilation scheme In general, an 01 analysis field, W(n) (denoting a vector at time t, with components c&(t,) where the subscripts i, j, k indicate the gridpoint (xi, yj, z,)), is the linear, least-squares best estimate of a prognostic field formed from a correction to the dynamical model forecast field @(n>:
explained
where K, =PH;(H,PH;+R)-I. H, defines the observational scheme at time t,; the superscript T denotes the matrix transpose; P is the forecast error covariance matrix and R is the observational error covariance. The error covariance of the analysis field, P,“, is P,” = (I-K,H,)P,
mode 0
Fig. 3. Statistics of the model ocean and the reference ocean for each mode. For the mean modal amplitude contours, the dashed line indicates negative values and the contour interval is 5 x 10’ m ’ s -’ for all modes. The contour interval for the modal amplitude standard deviation is 2.5 X lo3 m’s_’ for modes 0 and 1, and lo3 m2sC1 for mode 2.
M.M. Rienecker, D. Adamec /Journal
0.’
”
v
’
1000 -
x JI Ii $
2000.’ 3000-
*’
,’
.=’
3 4000-
5000 Fig. 4. The vertical eofs of streamfunction mined from the model ocean.
variability
deter-
where I is the identity matrix. The forecast error and the observational error are assumed to be uncorrelated. Details of this procedure are given by Miller (1986). In this study, the 01 was performed on a 2-dimensional field, so that au(n) is comprised of +t(t,) (elements on a 2-dimensional horizontal grid), and the covariance matrices are 2-dimensional. Two different methods to assimilate the surface altimeter height are considered here. For each, the upper layer streamfunction is used as a proxy for the surface altimeter height, assuming the relation: h altimeter
=
**f/g.
Throughout, the same assumptions were made with respect to the error covariances, viz., the measurement error was assumed to be white there is no correlation between measurement errors at different locations and no correlation between measurement errors and forecast errors. In addition, the model forecast error covariance, P was chosen to have an isotropic, homogeneous functional form similar to the observed streamfunction covariance:
P(r) = d[ 1 - (+)“I with a = 140 km and
e-(r’@, b =
(4)
120 km. The decorrela-
of Marine Systems 6 (1995) 125-143
131
tion scales used for this covariance are commensurate with those of the model streamfunction anomaly and assumes that the dominant contribution to the error in the forecast streamfunction comes from inadequate knowledge of the intensity and phase of the mesoscale variability. Phillips (1986) points out that this would be the case for an excellent forecast model that assimilates very little data. Certainly the cross-track updates fit into this scenario since (except in the region of the mid-latitude jet) the baroclinic information propagates across only a few grid points between updates and the update interval for a single track is about the same as the decorrelation time for the mesoscale eddy field (20-30 days). The assumed error covariance correlation scales are also commensurate with the cross-track spacing of the satellite orbits, rather than the along-track observation spacing, and is consistent with the findings of White et al. (1990b) and Rienecker and Miller (1990, i.e., if an isotropic error covariance function is used, the relevant decorrelation scale is no smaller than the larger of the meridional or zonal observation spacing. The reason for this, as Phillips (1986) points out, is that the error covariante scales estimated from observations are limited by the distance between the observations there is no information on finer scales. Only observations within a radius of 120 km were allowed to influence the analysis at any grid point since streamfunction data separated by 120 km have a correlation of only 0.1. This localization, the use of a single type of observation and the assumption of a spatially homogeneous observation noise level at a specified fraction of the model error variance mean that u2 can be an arbitrary normalization factor. It would be trivial to include a temporal variation for u2 which includes an assumed model error growth rate and the effect of reduced model forecast error variance at the points which have been corrected recently through assimilation (e.g., Rienecker and Miller, 1991). However, this has not been done here. The observational strategy is such that H,, is a matrix with non-zero (unit) entries only at model grid points which lie beneath a satellite orbit (see White et al., 1990a for a typical sampling strategy on the model grid).
132
M.M. Rienecker, D. Adamec /Journal
It should be pointed out that the use of approximate error covariance functions reduces the method to being a statistical interpolation rather than a statistically optimal interpolation. Daley (1991) discusses the impact of errors in assumed covariance functions. In practice, the use of assumed funtional forms for error covariances is usually the best that can be done and is still found to be effective. The two different techniques considered here assimilate: (1) the estimated amplitude, (Y~(x~,Yj, t,), of the dominant vertical eof of streamfunction variability, so that ~ij(t,> = (~i(x~, Yj, t,); and 2) the surface streamfunction, so that ~$~~(t,)= @i(xj, Yjy t,). Th e f irs t is referred to hereafter as the eof-based method, the second as the surfacebased method.
of Marine Systems 6 (1995) 125-143
The eof-based method
The assimilation cycle estimates (bG(t,) = c#xj, .Yj, f,J for each observation by assuming that all of the observation signal is contained in the first eof: (y;)(xi? Yj? ‘n) = Ilrf(‘i? Yj? ‘n)/511,
(5)
where tlk is the k-th component of the I-th eof of streamfunction variability and (cl$xi, yj, t,) is the surface layer streamfunction (surface height) observation sampled from the reference ocean. Since the eofs form a complete set of basis functions for the 3-layer streamfunction, the amplitude of the first eof from the model [c##,)] is obtained by projection: iv=3 af(XiY
Yj7 ‘n) =
C
(clk(xi,
k=l
Yj, fn)51k.
(6)
DAY 5000 Reference
Difference
layer 1
layer
2
layer 3
Fig. 5. The initial assimilation streamfunction analyses, the streamfunction from the reference ocean and the difference between the two, for each layer from one of the experiments conducted. The contour interval is 2 X lot, lo4 and lo4 m*s-’ in layers 1, 2 and 3, respectively. Contours less (greater) than - (+) 4 X lo4 in layer 1, and - (+) 2 X lo4 in layers 2 and 3 are shaded blue (orange).
M.M. Rienecker, D. Adamec /Journal
two-gyre experiment, but also marked differences in the meandering pattern of the mid-latitude jet and the intensity and placement of the eddies. The differences extend throughout the water column with similar scales (up to about 700 km). The inital root mean square difference (error> between the analysis and the reference field is 64% of the streamfunction variance of the reference ocean in the upper layer, 108% in the second layer and 89% in the bottom layer (the nrms statistic in Table 3). The seemingly high correlation (0.81 in the surface layer, 0.45 in the second layer and 0.58 in the bottom layer) show the similarities expected in the streamfunction distribution associated with the time-invariant surface forcing specified (only the basin-wide mean has been removed), as well as the influence of the bottom topography in steering the flow. After a complete orbit cycle (day 5020, Figs. 6a, 7a), both the eof method and the surface method have produced analyses much closer to the reference ocean, with differences smaller in scale (now about 200 km) throughout the water column for the eof method. The differences have been noticeably reduced in the far field of the mid-latitude jet (as also found by White et al., 1990a). The global rms differences have decreased by 52%, 41% and 11% of their initial values in each of the three layers, respectively, for the eof-based analyses (Table 3). In contrast, the rms error has decreased by 38% for the surface layer, but increased slightly in the subsurface layers for the surface-based analyses, and the differences are slightly larger scale. After three
After the 01 analysis is produced, the model integration proceeds with I,!J~
N=3 +
C af(Xi, 1=2
Yj9 fn)5[k,
133
of Marine Systems 6 (1995) 125-143
(7)
with (Y{,1= 2,3, determined analogous to Eq. (6). A variation on this would partition the observation signal across all eofs according to the variance explained by each, as in Dombrowsky and De Mey (1992); however since the higher eofs explain less than the assumed noise level (10% of the signal variance), such an approach is not warranted here. The inclusion of higher eofs would undoubtedly require a different form, or different parameters, for P in Eq. (4). The surface-based method
For this method, the 01 operates directly on the surface streamfuntion, $i(x,, yj, t,>. After the 01, the model integration proceeds with the analysis streamfunction in the surface layer, but the subsurface streamfunction unchanged.
3. Results
A comparison of the initial field used for one of the assimilation experiments (beginning model day 5000, Fig. 5) with the corresponding reference field from the series used for simulated data shows the similarities expected of this type of
Table 3 Root mean square error and correlation measures for one experiment. The root mean square difference between the analyses and the reference ocean is normalized by the global variance in the reference ocean (nrms) and given as a percentage Day
0 20 50 100
nrms C nrms C nrms C nrms C
EOF method
Surface method
Layer 1
Layer 2
Layer 3
Layer 1
Layer 2
Layer 3
64 0.81 31 0.96 32 0.96 28 0.97
108 0.45 64 0.84 71 0.81 51 0.87
89 0.58 75 0.75 80 0.76 55 0.86
64 0.81 40 0.92 44 0.90 58 0.83
108 0.45 115 0.48 141 0.54 239 0.26
89 0.58 101 0.53 132 0.62 250 0.45
MM. Rienecker, D. Adamec /Journal
134
of Marine Systems 6 (1995) 125-143
DAY 5020 Difference
layer
1
layer
2
layer 3
DAY 5050 Difference
layer 1
layer 2
layer 3
(a)
M.M. Rienecker, D. Adamec /Journal
of Marine Systems 6 (1995) 125-143
135
Difference
layer 1
layer 2
layer 3
Fig. 6. The analyses from the eof-based assimilation technique, streamfunction from the reference ocean and the difference between the two, for each layer. The fields after (a) 20 days (b) 50 days and (c) 100 days of assimilation. Contour intervals and shading as for Fig. 5.
repeat cycles (day 5050), it appears that the similarity between the eof-based analysis and the reference ocean has not markedly improved, except that the region of greatest differences is more tightly bound to the jet (Fig. 6b). In contrast, the differences in the surface-based analysis have increased in scale to about 1000 km in the subsurface layers and the magnitude of the differences has increased substantially, to (globally) 141% of the natural variability for layer 2 and 132% for layer 3 (Fig. 7b, Table 3). After 100 days of assimilation, the eof-based analyses are visually in very close agreement with the reference ocean in all layers (Fig. 6c) and the global statistics confirm the good agreement in all layers (Table 3), with a reduction of the initial error by 56% in the surface layer, 53% in the second layer
and 38% in the bottom layer. Whereas there are some features in the surface-based analysis which compare well with the reference field, there are also large discrepancies (Fig. 7c, note the difference in contour interval) and the correlation between the analyses and the reference field has decreased significantly (Table 3). In the surface layer, the intense discrepancies are still small scale, but in the subsurface layers the differences are larger scale (up to 1000 km) and are not confined to the region of the jet. This is evident in the comparisons of Fig. 8 where the absolute error of the analysis field is normalized by the standard deviation of the variability at each grid point in the reference ocean. After 100 days of assimilation, the improvement in the eof-based analysis is global, whereas the disagreement (note
M.M. Rienecker, D. Adamec/Journal
136
of Marine Systems 6 (1995) 125-143
DAY 5020
layer
1
layer
2
layer
3
DAY 5050 Reference
layer 1
layer
2
layer
3
Difference
M. M. Rienecker, D. Adamec /Journal
of Marine Systems 6 (1995) 125-143
137
DAY 5100 Difference
layer 1
layer 2
layer 3
Fig. 7. The analyses from the surface-based assimilation technique, streamfunction from the reference ocean and the difference between the two, for each layer. The fields after (a) 20 days (b) 50 days and (c) 100 days of assimilation. The contour intervals are as for Fig. 5, except for 100 days when the contour intetval in layers 2 and 3 are the same as for layer 1.
again the difference in contour interval for the surface- based method at day 5100) in the subsurface layers of the surface-based method extend throughout a large fraction of the domain, but are most intense along the western boundary. The large relative amplitude along the eastern boundary is merely due to the low variance (normalizing factor) there, as seen in Fig. 3. In this region. linear dynamics apply and the application of a correction to the forecast streamfunction is akin to wind stress curl forcing of Rossby waves which are then clearly seen in the surface signal. The time series of the root mean square difference between the analyses and the reference ocean for these experiments are shown in Fig. 9a. Most of the reduction in error occurs very quickly, over the first satellite repeat cycle. After two
repeat cycles (34 days), the error level for the eof-based analyses is fairly constant, with excursions about the asymptote due to mesoscale activity not resolved by the satellite observations. Finer observation spacing would lower the error bound, as shown by Holland and Malanotte-Rizzoli (1989). The error level in the bottom layer slowly decreases in time, and it would likely take much longer (gyre recirculation timescales) for the error level to reach its minimum. In contrast, the error levels for the surface-based method decrease dramatically only in the surface layer. The errors stay low for about 50 days in all layers, but then start to increase without any apparent bound. Eventually, as found by White et al (1990a), the initial model field would have given a closer agreement with the reference ocean than the
MM
138
Rienecker, D. Adamec /Journal
qf Marine Systems 6 (1995) 125-143
Day 5020
EOF Method
Surface
Method
EOF Method
Day 5100 Surface
Method
layer 1
layer
2
layer
3
Fig. 8. The absolute difference between the analyses and the reference ocean, normalized by the standard deviation of the streamfunction variability at each point in the reference ocean. The contour interval is 2, except for the last column where the contour interval is 4. The shading is blue for values between 2 and 4 (4 and 8 for the last column) and orange for values greater than 4 (greater than 8 for the last column).
analysis field and the model should be re-initialized. The results for this single experiment is typical of that for each experiment in the ensemble (Fig. 9b), indicating the results are robust, at least for the dynamical scenario (wind forcing, bottom topography, stratification, resolution) considered here. The rrns error statistics for the dynamical modes (Fig. 912,d) explain the differences in behaviour of the two methods. The eof-based method shows a marked decrease in the error in the barotropic and first baroclinic modes and a slight decrease for the second baroclinic mode
over the entire assimilation experiment. In contrast, the surface-based method shows initial improvement in the barotropic and first baroclinic modes, but subsequently, the error is reduced only in the first baroclinic mode for which the error remains low for about 100 days. The error in the barotropic mode increases markedly after one or two satellite repeat cycles and the error in the second baroclinic mode increases throughout the entire assimilation experiment. The explanation for the differences in modal errors lies in the decomposition of the correction, ef(n), to the model forecast streamfunction. For
M.M. Rienecker, D. Adamec /Journal
of Marine Systems 6 (1995) 125-143
139 (b)
.-’ .A____ r”....“‘.-...........~..................,...................~. _____________________
5.0-103-~L.,..
50
100 assimilation day
0: 0
150
I 50
I 100 assimilation day
7
150
k.4
ot. 0
1.
50
I”‘.1
,‘,
150
100 assimilation day
ot 0
.
I
.
.
.
.
50
.
100
.
...1 150
assimilation day
(4
:..................................................,,..................... ===z
~____-_---r=========;; !
0
50
,...
I
100 assimilation day
I.
150
Fig. 9. Time series of root mean square error during the assimilation cycle. (a) A single assimilation experiment, the solid line is the streamfunction error in layer 1, the dashed line is the error in layer 2 and the dotted line is the error in layer 3. The thicker lines are for the eof-based method, the thinner lines for the surface method. (b) As for (a), but the ensemble average error. (c) As for (a), but for the modes. The solid line is for the barotropic mode, the dashed line for the first baroclinic mode, and the dashed line for the second baroclinic mode. Cd) As for Cc) but for the ensemble average error. (e) As for (a) but for potential vorticity (scaled by 1O6).
M.M. Rienecker, D. Adamec /Journal
140
the eof-based method, 28.2% of the variance goes to the barotropic mode, 64.8% to the first baroclinic mode and 7.0% to the second baroclinic mode, very close to the distribution for the reference ocean. In contrast, for the surface-based method, only 6% of the variance goes to the barotropic mode, 65% to the first baroclinic mode and 29% to the second baroclinic mode. The continuous addition of anomalous energy in the higher baroclinic mode apparently leads to larger transfers to barotropic energy through the various energy cascades, some of which may be enhanced by topography (e.g. Adamec, 1988). The larger subsurface scales in the difference fields for the surface-based method (Fig. 7) are associated with the error in the barotropic mode and its scales as set by topography. For the surface layer, the larger scales are well-defined by the observation strategy (viz., scales of about 300 km for the Geosat track spacing, as shown by White et al., 1990b), so the errors remain in the smaller (unresolved) scales. It is of interest to compare the present vertical projection technique with that used by Haines (1991) wherein the streamfunction in the surface layer is replaced with the observed streamfunction and the subsurface streamfunction is adjusted so that the subsurface potential vorticity, &, is unchanged, where
fll
!L= v*V+,+f(y)
- #Q-i,*
-&Cc,,*)
k +s
f&B i’.‘,k HN .
This process entails solving elliptic equations for each subsurface streamfunction correction: V26$, - (Y& + Y22,,)& + Y&V3 = -Y;,,~&~ V2S6, -
Y32,2W3
+
d,*w2
=
0,
where, using Haines’ notation, yi k +, = ft/g; +1,2Hk and c?$~ is the streamfunction correction in layer k. With a spectral (horizontal) decomposition and the model parameters used in this study, it is easy to show that for very small wavenumbers, the response is S& = ?j+2 = S+i, i.e., the constant potential vorticity method projects the surface signal onto the
of Marine Systems 6 (1995) 125-143
barotropic mode. For wavenumbers = (70 km)-‘, the response is 8$s = 0.2564, and SI&*2: 0.46$,, ratios which are very similar to those for the first eof in this study which has 5,s = 0.205,, and 5i2 = 0.375ii. For smaller scales, the surface signal has much reduced influence on the subsurface layers and at the first baroclinic radius of deformation, the influence is very small, particularly in the bottom layer: S$, = O.O6S$, and &,!I, = 0.2s*,. In an identical-twin, flat bottom ocean study, Haines found that the error in the bottom streamfunction increased during assimilation but that afterwards the bottom Ekman drag and the interfacial eddy stresses force the bottom layer towards the control ocean. The damping effect of bottom drag and the “correcting influence of the currents generated from above” help redistribute the potential vorticity in the bottom layer. However, the dynamical coupling is two-way; i.e., the variability (errors) in the bottom layer also force variability in the middle layer(s) (see Lozier and Riser, 1990, for the long-term balances in eddy driving of subsurface layers in a flat-bottom ocean) and the errors may magnify in time. In addition, because of the thickness of the bottom layer, errors in potential vorticity in this layer, either in the initial conditions or introduced through assimilation, have predominant influence on the barotropic mode. This mode has the strongest interaction with the bottom topography. Outside of smooth, large-scale bottom slopes, topographic interactions tend to change the energy cascade toward smaller scale barotropic flow as well as surface-intensified baroclinic flow (e.g. Adamec, 1988), the latter requiring further adjustments in the upper layer flow in the assimilation process. For the eof method, the surface information is projected onto the subsurface layers according to the long-term statistical relationship between the layers - these statistics automatically encompass the energy cascades consistent with the underlying dynamics, including the surface and bottom interactions, and any statistically significant potential vorticity homogenization. Of course, the streamfunction ratios are not correct at all places and time. For both the constant potential vortic-
M.M. Rienecker, D. Adamec /Journal
ity method and the eof method, the reliance is on the model to adjust the incorrect potential vorticity according to the new analysis streamfunction in each layer and the dynamical balances of the model. The difference between the potential vorticity for the analyses and the reference ocean is shown in Fig. 9e. For this case study, the potential vorticity difference in the lower layers remains fairly constant for the eof-based method. The assimilation procedure reduces the potential vorticity error in the upper layer by 39% after 150 days and increases the potential vorticity error by 15% in the middle layer and by 0.4% in the bottom layer. Maps of potential vorticity (not shown) show that the large pool of fairly homogeneous potential vorticity in the middle layer, beneath the surface mid-latitude jet, is not destroyed by the eof-based analyses. The bottom layer has no such homogeneous region because of the bottom topography.
4. Discussion It has been demonstrated that optimal interpolation and the use of a single model-based vertical eof of streamfunction variability are effective and efficient for assimilating surface altimeter data into fraternal-twin simulations of a 3-layer quasigeostrophic model with topography. After 100 days of assimilation (about 6 satellite repeat cycles), the error in the upper two layers has been reduced by over 50% and that in the bottom layer by 38%. The essence of the method is that the eofs capture the statistics of the dynamical balances in the model and ensure that this balance is not inappropriately disturbed during the assimilation process. The analysis fields have been filtered according to the model’s dynamical and statistical balances. In contrast, a method which assimilates the information at the surface only and relies on the model to extend the information throughout the water column imposes an incorrect distribution of energy on the fundamental dynamical modes of the model; the reliance then is on the model dynamics to redistribute the energy according to its own nonlinear cascade. Of course, if the assimilation were curtailed, the
of Marine Systems 6 (1995) 125-143
141
model would do exactly that, reproducing, after a few eddy turn-over times, an energy distribution more like that in the model ocean. This suggests that there may be some benefit (A. Robinson, pers. commun.) in assimilating data infrequently and allowing the model to redistribute the energy according to its own dynamical balances. For the eof-based method, such considerations are not necessary and the data may be used close to its observation time to try to preserve the phase information as closely as possible. The eof-based method effectively reduces the degrees of freedom and so the computational burden of the assimilation process. The success of using a single eof for the entire domain relies somewhat on the fact that here this single eof explains more than 90% of the entire lateral and vertical streamfunction variance. Such high levels of variance explained by a single mode are not unprecedented in data: observations in limited area domains (e.g., Rienecker et al., 1987; Dombrowsky and De Mey, 1992) show dominance by a single local vertical mode of dynamic topography. In practice, the implementation of additional and spatially varying eofs would pose no difficulties to the use of this technique for effective assimilation of surface information. The implementation requires the specification of a horizontal forecast error covariance function for each eof. The underlying assumption is that the 3-dimensional error covariante is separable (horizontal from vertical) so that the analysis may be decoupled in a similar manner to the vertical decoupling of the dynamical equations by projection onto dynamical vertical modes. The use of eofs (statistical modes) is more useful than dynamical modes if (as in this study) the variance is predominantly in a few statistical modes. For special forms of the vertical observation error covariance matrix (see Daley, 1991, Ch. 41, it can be shown that the eofs of the background error covariance (here assumed to be the same as the eofs of streamfunction variability) are also eofs (and so an efficient characterization) of the analysis error covariance. The full Kalman filter would evolve this analysis error covariance according to the nonlinear, vertically-coupled dynamics. However, as long as the system (model) error is not dominant, it might be expected that
142
M.M. Rienecker, D. Adamec /Journal
(unless a lot of data is assimilated often) the initial eofs would remain a good representation of the error covariance as they are of the streamfunction variability which has evolved under the same dynamics. Even with a model with much higher vertical resolution, it is unlikely that a large number of eofs would be required. The use of a small number of eofs of the forecast error covariance matrix has been found efficacious in the National Meteorological Center’s global data assimilation system (Parrish and Derber, 1992). For the ocean, in a joint eof decomposition of fields of temperature, salinity, oxygen and nutrient observations in the North Atlantic, Fukumori and Wunsch (1991) were able to explain 92% of the variance with 6 modes, 86% with 4 modes and 68% with only 2 modes. For model-based eofs, the explained variance for the low modes is likely to be higher than for observations because of lower noise levels. Lack of knowledge of surface forcing parameters and uncertainty in fundamental model parameters could lead to a model-estimated eofs which are substantially different from that supported by analysis fields influenced by a lot of data (because of different ratios of barotropic to baroclinic energy etc). For this case, new eofs could be reestimated after many assimilation cycles. In the assimilation of satellite observations (rather than simulated data) it will be important to include along-track observational error covariances.
Acknowledgements The authors gratefully acknowledge enlightening discussions with and suggestions from Mike Bell, Bob Haney, Robert Miller, Allan Robinson, and John Marshall. The comments of two anonymous reviewers helped clarify the presentation and are greatly appreciated. This research was supported by the Ocean Modeling and Data Analysis Program under NASA’s Climate and Hydrologic Systems Branch. Computing was carried out at the NASA Center for Computational Sciences at Goddard Space Flight Center.
of Marine Systems 6 (1995) 125-143
References Adamec, D., 1988. Numerical simulations of the effects of seamounts and vertical resolution on strong ocean flows. J. Phys. Oceanogr., 18: 258-269. Bengtsson, L., Ghil, M. and Kallen, E. (Editors), 1981. Dynamic Meteorology: Data Assimilation Methods. Springer, New York, 330 pp. Cummins, P. and Mysak, L.A., 1988. A Quasi-geostrophic circulation model of the Northeast Pacific, Part I: A preliminary numerical experiment, J. Phys. Oceanogr., 18: 1261-1286. Daley, R., 1991. Atmospheric Data Analysis. Cambridge University Press, 457 pp. De Mey, P. and Robinson, A.R., 1987. Assimilation of altimeter eddy fields in a limited-area quasigeostrophic model. J. Phys. Oceanogr., 17: 2280-2293. Dombrowsky, E. and De Mey, P., 1992. Continuous assimilation in an open domain of the Northeast Atlantic 1. Methodology and application to AthenaA-88. J. Geophys. Res., 97: 9719-9731. Fukumori, I. and Wunsch, C., 1991. Efficient representation of the North Atlantic hydrographic and chemical distributions. Prog. Oceanogr., 27: 111-195. Haines, K., 1991. A direct method for assimilating sea surface height data into ocean models with adjustments to the deep circulation, J. Phys. Oceanogr., 21: 843-868. Holland, W.R., 1989. Altimeter data assimilation into ocean circulation models - some preliminary results. In: D.L.T. Anderson and J. Willebrand (Editors), Oceanic Circulation Models: Combining Data and Dynamics. Kluwer, Dordrechts, 605 pp. Holland, W.R. and Malanotte-Rizzoli, P., 1989. Assimilation of altimeter data into an ocean circulation model: space versus time resolution studies. J. Phys. Oceanogr., 19: 1507-1534. Hurlburt, H., 1986. Dynamic transfer of simulated altimeter data into subsurface information by a numerical ocean model. J. Geophys. Res., 91: 2372-2400. Lozier, MS. and Riser, S., 1990. Potential vorticity sources and sinks in a quasi-geostrophic ocean: beyond Western Boundary Currents. J. Phys. Oceanogr., 20: 1608-1627. Marshall, J. and Berry, P., 1991. Ocean modelling studies inspired by the Seasat altimeter. Int. J. Remote Sensing, 12: 1631-1638. Mellor, G.L. and Ezer, T., 1991. A Gulf Stream model and an Altimetry assimilation scheme. J. Geophys. Res., 96: 8779-8795. Miller, R.N., 1986. Toward the application of the Kalman filter to regional open ocean modeling. J. Phys. Oceanogr., 16: 72-86. Parrish, D.F. and Derber, J.C, 1992. The National Meteorological Center’s Spectral Statistical-Interpolation Analysis Scheme. Mon. Weather Rev., 120: 1747-1763. Phillips, N.A., 1986. The spatial statistics of random geostrophic modes and first-guess errors. Tellus, 38A: 314-332.
M.M. Rienecker, D. Adamec /Journal
Rienecker, M.M. and Miller, R.N., 1991. Ocean data assimilation using Optimal Interpolation with a quasi-geostrophic model. J. Geophys. Res., 96, l&093-15,103. Rienecker, M.M., Mooers, C.N.K. and Robinson, A.R., 1987. Dynamical interpolation and forecast of the evolution of mesoscale features off Northern California. J. Phys. Oceanogr., 17: 1189-1213. Verron, J., 1990. Altimeter data assimilation into an ocean model: sensitivity to orbital parameters. J. Geophys. Res., 95, 11,443-11,459. Verron, J., 1992. Nudging satellite altimeter data into quasi-
of Marine Sjstems 6 (1995) 125-143
143
geostrophic ocean models. J. Geophys. Res., 97: 74797491. White, W.B., Tai, C.-K. and Holland, W.R., 1990a. Continuous assimilation of simulated Geosat altimetric sea level into an eddy-resolving numerical ocean model, 1, Sea level differences. J. Geophys. Res., 95: 3219-3234. White, W.B., Tai, C.-K. and Holland, W.R., 1990b. Ccmtinuous assimilation of simulated Geosat altimetric sea level into an eddy-resolving numerical ocean model, 2, Referenced Sea level differences. J. Geophys. Res., 95: 32353251.