A simple method for reducing seasonal bias and drift in eddy resolving ocean models

Ocean Modelling 13 (2006) 109–125 www.elsevier.com/locate/ocemod

A simple method for reducing seasonal bias and drift in eddy resolving ocean models Keith R. Thompson a

a,*

, Daniel G. Wright b, Youyu Lu b, Entcho Demirov

c

Department of Oceanography, Dalhousie University, Halifax, NS, Canada b Bedford Institute of Oceanography, Dartmouth, NS, Canada c Department of Physics, Memorial University, St. Johns, Nfld, Canada

Received 29 July 2005; received in revised form 6 October 2005; accepted 30 November 2005 Available online 9 January 2006

Abstract An effective and computationally efficient method is described for reducing bias and drift in eddy resolving ocean models. The basic idea is to nudge the model towards gridded climatologies of observed temperature and salinity in prescribed frequency-wavenumber bands; outside of these bands the model’s dynamics are not directly affected by the nudging and the model state can evolve prognostically. Given the restriction of the nudging to certain frequency-wavenumber bands, the method is termed spectral nudging. The frequency-wavenumber bands are chosen to capture the information in the climatology and thus are centered on the climatological frequencies of 0, 1 cycle per year and its harmonics, and also low wavenumbers (consistent with the smooth nature of gridded climatologies). The method is first tested using a linear, barotropic ocean model forced by a time-varying wind stress curl. For this example it is possible to give explicit expressions for the effect of spectral nudging on the model’s response to wind forcing. The method is then applied to an eddy permitting model of the North Atlantic driven by realistic surface fluxes. It is shown that the model maintains a statistical steady state over the several decades of integration with no evidence of bias or drift. Further the seasonal cycle of coastal sea level, and the spatial distribution of sea level variance, are shown to agree well with independent sea level measurements made by tide gauges and altimeters, respectively. The application of spectral nudging to shelf and coupled atmosphere–ocean modelling is discussed.  2005 Elsevier Ltd. All rights reserved. Keywords: Ocean models; Bias; Drift; North Atlantic; Sea level

1. Introduction Bias and drift are common problems with ocean models (e.g. Bleck, 2002; Halliwell, 2004; Huddleston et al., 2004). Many factors can contribute to bias and drift including inadequate model resolution, poor parameterizations of sub-grid scale processes and inaccurate surface and lateral boundary conditions. Given

*

Corresponding author. E-mail address: [email protected] (K.R. Thompson).

1463-5003/$ - see front matter  2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.ocemod.2005.11.003

110

K.R. Thompson et al. / Ocean Modelling 13 (2006) 109–125

the range and complexity of potential contributors, it is hardly surprising that a general solution to these problems is not available. The objective of this study is to develop a simple method for suppressing bias and drift in eddy resolving models when the variability of interest has periods of days to seasons. The method proposed here may be used in the reconstruction of historical changes of the ocean state and also making forecasts with lead times of days to seasons. The method is based on the simple idea of nudging the model toward seasonal climatologies of ocean observations in specified frequency and wavenumber bands; outside of these bands the model is unconstrained. The Bayesian framework (e.g. Lorenc, 1986) is arguably the most useful in terms of motivating the new method, and relating it to other assimilation approaches. Let xt+1 = /t(xt, wt) represent a nonlinear dynamic model that describes how the discretized ocean state vector, xt, is advanced forward one time step. In one of our applications the nonlinear dynamic operator /t corresponds to a fully nonlinear, eddy permitting model of the North Atlantic. Uncertainty in the forecast is reflected by the inclusion of the model error vector wt. The relationship between xt and the contemporaneous vector of observations, yt, is modelled by yt = ht(xt, vt) where ht is the observation operator and vt defines the observation error (e.g. Daley, 1994). It is well known (e.g. Jazwinski, 1970) that if the error terms are mutually independent, and Yt denotes all available observations up to and including time t, then routine application of Bayes’ R Rule allows the probability distribution of xt given Yt to be updated sequentially using pðxt jY t Þ / pðy t jxt Þ pðxt jxt1 Þpðxt1 jY t1 Þ dxt1 . This conditional density contains all the information on xt given Yt; quantities like the conditional mean of xt, its covariance and higher order moments can all be obtained from this probability density function. Unfortunately the updating equation given in the previous paragraph can only be solved by direct integration for low dimensional problems (e.g. Kitagawa, 1998); approximate numerical schemes must be developed for use with high-dimensional ocean models. If /t and ht are linear operators, and wt and vt are normally distributed with zero mean, the celebrated linear Kalman Filter (e.g. Jazwinski, 1970) provides an elegant way of sequentially updating the conditional mean and variance of xt given Yt. For nonlinear problems the situation is less clear cut and a number of sub-optimal schemes have been proposed including the extended Kalman Filter, the singular evolutive extended Kalman Filter (Pham et al., 1998) and the ensemble Kalman Filter (Evensen, 1994). All forms of Kalman Filter correspond to ‘nudging’ the ocean state toward the observations using an equation of the form xt ¼ xft þ K t ðy t  ht ðxft ÞÞ where xft ¼ /t ðxt Þ is the one-step-ahead-forecast, and Kt are the nudging coefficients organized in the so-called Kalman gain matrix. The differences in the various Kalman Filters affect only the calculation of Kt. A major simplification is to prescribe Kt based on an assumed covariance structure of the forecast errors (e.g. Daley, 1994). If Kt is assumed to be diagonal then the nudging terms take the conventional form cðy t  xft Þ at observation points where c1 corresponds to a relaxation time in model time steps. Sarmiento and Bryan (1982) used this form of nudging to assimilate observations of temperature and salinity into a realistic ocean model. Their technique has proved very useful over the years and is usually referred to as the robust diagnostic method. If c = 1 the nudged ocean state is replaced by directly observed values, a procedure known as direct insertion. It corresponds to the standard diagnostic method of Sarkisyan (1977) if temperature and salinity data are assimilated. The main advantages of the standard and robust diagnostic forms of nudging are that they are easy to implement, robust and can keep the model arbitrarily close to the observations. Unfortunately there are serious limitations including the suppression of eddies when nudging toward climatological observations, and the introduction of artificial phase lags in the model response (e.g. Woodgate and Killworth, 1997). An important assumption underlying the Kalman Filter is that the model’s predictions are unbiased. Direct comparison of time-averaged model predictions and observations suggests that this is often not the case. It is true that sequential schemes like the Kalman Filter, which are based on estimates of forecast error variability about the mean, can still be used to nudge a biased model toward observations but this will be achieved in a suboptimal way; it is highly unlikely that the covariance structure of the forecast errors will match the spatial structure of the bias errors. (Thacker and Esenkov (2002) provide a good discussion of this point.) Our solution to this problem is to suppress the bias before the observations are assimilated. In this paper we outline a straightforward method that ensures the model’s climatology does not stray too far from an observed climatology, thereby suppressing bias and drift. The basic idea is to nudge the model in frequency and wavenumber bands that cover the frequencies and wavenumbers that are well-resolved in the

K.R. Thompson et al. / Ocean Modelling 13 (2006) 109–125

111

climatology. In accord with the naming of a related technique used to downscale large-scale atmospheric states using regional models (von Storch et al., 2000 and references therin) we refer to the technique as spectral nudging. Outside of these bands the nudging is effectively zero and the model can evolve prognostically. This contrasts sharply to more conventional nudging schemes which suppress eddy variability more strongly as the strength of the nudging increases (e.g. Woodgate and Killworth, 1997). Spectral nudging is designed to be implemented sequentially and so it is possible to assimilate conventional data, such as in situ temperature and salinity or sea level from altimeters, in the usual way (e.g. using Kalman Filters or the approach of Cooper and Haines (1996) for altimeter data) after the suppression of the bias and drift by the spectral nudging step. To illustrate spectral nudging we introduce in Section 2 a highly idealized, linear ocean circulation model forced by a time-varying wind stress curl. In Section 3 we give explicit expressions for the effect of conventional nudging on the dynamics of this linear model. Spectral nudging of the model is fully described Section 4, including explicit expressions for its effect on the frequency-dependent transfer functions linking the wind forcing and ocean response. The application of spectral nudging to a 1/3 degree model of the North Atlantic is described in Section 5 and results are compared against sea level observations made by coastal tide gauges and altimeters. Results are summarized in the final section where suggestions are made for future work, including spectral nudging of shelf and coupled atmosphere–ocean models. 2. Idealized ocean model We now introduce a highly idealized, linear barotropic ocean model in order to illustrate the limitations of conventional nudging and also to indicate a simple way of overcoming them. One advantage of using such a simple model is that it is possible to obtain explicit expressions for the effect of nudging on the model’s dynamics and the ocean response. The application of spectral nudging to fully nonlinear ocean models, where variations in different frequency-wavenumber bands can interact, is discussed in Section 5. The idealized ocean model has a rigid lid, a flat bottom and lies on a mid-latitude b-plane. Its stream function w evolves according to the following forced Rossby wave equation (e.g. Gill, 1982):   o o2 w ow o2 w ¼ d þF ð1Þ þ ot ox2 ox ox2 where x is the east–west coordinate, non-dimensionalized by the ocean width L, and t is time, non-dimensionalized by (Lb)1 which is order 1 day for a typical mid-latitude basin. The parameter d is the width of the western boundary current region scaled by L (i.e. the non-dimensional Stommel width); its reciprocal is the spin-down time of the system in non-dimensional time. F denotes forcing by the wind stress curl and is allowed to vary with both x and t; its characteristic magnitude sets the scale for w. The coastal boundary conditions are w = 0 along the west (x = 0) and east (x = 1) boundaries of the model domain. The discretized form of (1) can be written in the form wt ¼ Uwt1 þ Qat

ð2Þ

where wt is the n · 1 state vector at time t, U is the n · n dynamics matrix that depends on d and the space and time steps, and Qat is the n · 1 forcing vector. The m · 1 vector at determines the time-varying amplitude of m forcing components; the spatial structures of the forcing components are given by the columns of the n · m matrix, Q. Typical forcing and the stream function response are shown in Fig. 1. The leftmost panel shows the time-varying amplitude of the forcing components (i.e. the two components of at). One component (red line) is designed to generate a response at ‘‘eddy’’ time scales. It is a realization from an AR2 process (e.g. Priestley, 1981). The other component (black line) represents climatological forcing. The spatial structures of both forcing components (i.e. the two columns of Q) are identical and represent spatially uniform forcing in (1). The remaining three panels show the stream function response to the total forcing (Qat) and the two components separately. The linearity of the model implies that the second panel (total w) is the sum of the third and fourth panels. The response to climatological forcing is close to steady state, as expected for a barotropic model (e.g. Gill, 1982); it includes a western boundary flow that varies almost sychronously with the seasonal wind stress curl, in agreement with simple Sverdrup theory and the associated return flow. The eddy response is straightforward to interpret in terms

112

K.R. Thompson et al. / Ocean Modelling 13 (2006) 109–125

Fig. 1. Response of the idealized ocean model to stochastic and climatological forcing. The left panel shows the time variation of the forcing components (at). The remaining three panels show the response of the stream function, as a function of location and time, to the total forcing and individual forcing components. The climatological forcing is a unit amplitude sinusoid with a period of 1 year; the stochastic component generates variability at eddy time scales and is a realization from an AR2 process with a damped resonance centered on a period of 32 days. (The autoregressive coefficients are a1 =  1.9 and a2 = 0.94.) The spatial structure of both forcing terms is uniform. The friction parameter d was 0.04, corresponding to a spin down time of 25 days. The nondimensional model time step is 1, corresponding to 1 day in dimensional terms. All four panels are for the last year of a 10 year integration of the model. There are 200 active w grid points in the x-direction.

of westward propagating barotropic Rossby waves that are dissipated in the western boundary current region after reflecting as slowly moving, short Rossby waves (Pedlosky, 1965). To describe the model response in the frequency domain, assume the forcing at and streamfunction wt are in a statistical steady state and can be described using the following spectral representation: Z p at ¼ eixt daðxÞ ð3Þ p Z p wt ¼ eixt dwðxÞ ð4Þ p

where, loosely speaking, da(x) and dw(x) are vectors of complex amplitudes of the sinusoidal components of at and wt. (See Priestley (1981) for a complete discussion of the spectral representation of multivariate processes.) Substituting these expressions into (2) leads to the following relationship between forcing and response at frequency x: dwðxÞ ¼ T u daðxÞ

ð5Þ

where 1

T u ðxÞ ¼ ½I  Ueix  Q

ð6Þ

K.R. Thompson et al. / Ocean Modelling 13 (2006) 109–125

113

Fig. 2. Effect of nudging on the frequency response of the idealized model. For all panels the x-axis shows the east–west location and the y-axis shows the frequency of the forcing (in cycles per day). The upper panels show the gain (jTj, where T is the transfer function); the bottom panels show the phase (arg T, in degrees). The first column of panels is for no nudging (see (6)), the middle column is for conventional nudging ((9) with c = 0.1), and the last column is for spectral nudging ((14) with c = 1 at 0 and 1 cycle per year, j1 = 1000 d). Other model settings are identical to those defined in Fig. 1.

is the n · m matrix of transfer functions that relates the amplitude and phase of the m forcing components to the stream function response at the n grid points. The left column of panels in Fig. 2 shows the gain (jTuj) and phase (arg Tu) as a function of x and x for spatially uniform forcing. At low frequency (i.e. periods exceeding about 200 days) the response is amplified in the western boundary current region and approximately out of phase with the wind forcing, consistent with theory for simple, barotropic wind-driven circulation. At forcing periods between about 40 and 80 days, there is an amplified response that may be interpreted in terms of resonantly forced, basin modes (Longuet-Higgins, 1964). 3. Conventional nudging Let wct denote the observed stream function climatology with which the idealized ocean model is required to agree. In general wct will not coincide with the climatological model response to forcing by at. The simplest approach is to nudge the model as follows: wt ¼ wft þ cðwct  wft Þ where

wft

ð7Þ

is the one-step-ahead forecast:

wft ¼ Uwt1 þ Qat The last term in (7) nudges the ocean forecast toward the climatology with a relaxation time of c steps.

ð8Þ 1

model time

114

K.R. Thompson et al. / Ocean Modelling 13 (2006) 109–125

Fig. 3 shows the effect of nudging the model towards wct in the time domain. As expected, increasing the strength of the nudging term (i.e. increasing c) moves the solution from the prognostic solution (c = 0, left panel) toward the observed climatology (c = 1, right panel). As expected the nudging also suppresses the non-climatological variability. This is clearly evident in the following expression for the transfer function of the nudged system for frequencies at which the climatology has no energy (i.e. excluding 0 and 1 cycle per year in this case): 1

T n ¼ ð1  cÞ½I  ð1  cÞUeix  Q

ð9Þ

Note that if c = 0 in this equation then Tn = Tu; if c = 1 the response is zero. At intermediate values of c, it is clear that the nudging introduces artificial phase lags and amplitude changes that distort the dynamics. This is illustrated in Fig. 2 which shows the gain and phase of the nudged system. It is clear from this figure that the nudging has not only suppressed the dominant basin mode, it has also changed the westward propagation speed of information across the basin. This is consistent with an earlier study by Woodgate and Killworth (1997) that showed nudging the shallow water equations can seriously alter the propagation of inertia-gravity waves. One way of removing the artificial phase lags in this linear system is to replace (7) by wt ¼ wft þ cðwct  wft  Uðwct1  wt1 ÞÞ

ð10Þ

In this formulation the nudge is proportional to the difference between the climatology and model forecast at time t (i.e. wct  wft ) and the difference between the climatology and the analyzed state at the previous time step (i.e. wct1  wt1 ). The advantage of this new formulation is that the transfer function corresponding to (10) takes on a particularly simple form:

Fig. 3. Suppression of eddies by conventional nudging. Each panel shows the response of the stream function to the forcing defined in Fig. 1. Eq. (7) was used with c increasing from 0 (leftmost panel) to 1 (rightmost panel). The leftmost panel corresponds to no nudging and is identical to the second panel of Fig. 1; the rightmost panel corresponds to direct insertion of the observed climatology to which the model is nudged, i.e., wct ¼ sinðpxÞ sinð2pt=365Þ. Other model settings are identical to those defined in Fig. 1.

K.R. Thompson et al. / Ocean Modelling 13 (2006) 109–125

T n ¼ ð1  cÞT u ðxÞ

115

ð11Þ

where Tu(x) is the un-nudged transfer function given by (6). It is clear that the new formulation introduces no artificial phase lags into the dynamics; it simply scales Tu(x) by the constant factor 1  c. It is also straightforward to show that the (spun up) response in the time domain can be written wt ¼ ð1  cÞwut þ cwct

ð12Þ

where wut is the un-nudged solution of (2). Thus as c increases from 0 to 1, the stream function changes smoothly, and with no artificial phase lags, from the un-nudged solution (c = 0) to the observed climatology (c = 1). But the problem of attenuating eddy variability by the nudging remains. In the next section we use (11) to motivate a straightforward solution to this problem. To physically interpret (10) we note that, on using (8), it can be rewritten in the mathematically equivalent form wt ¼ wft þ cðwct  Uwct1  Qat Þ

ð13Þ

Thus for this linear error-free model, the nudge in (10) is proportional to the difference between the forcing supplied by the model (Qat) and the forcing required by the model to make the climatology evolve from its observed value at time t  1 to its observed value at time t (i.e. wct  Uwct1 ). 4. Spectral nudging One way of overcoming the suppression of eddy variability by conventional nudging is by making c in (11) a function of frequency: T s ¼ ½1  cCðxÞT u ðxÞ

ð14Þ

where Ts is the transfer function of the spectrally nudged model, and C(x) is a function of frequency yet to be specified. From (14) it is straightforward to obtain the following frequency-dependent generalization of (12): dwðxÞ ¼ ½1  cCðxÞ dwu ðxÞ þ cCðxÞ dwc ðxÞ

ð15Þ

From these two equations it is clear that if we require C = 0 at non-climatological frequencies, then Ts = Tu and dw = dwu; the nudging will have no effect on the dynamics or the ocean response, as required. If we require C = 1 at climatological frequencies (henceforth xc) then dw(xc) will be a simple linear combination of the un-nudged response, dwu(xc), and the observed climatology, dwc(xc). Further, if c = 1, the model will agree exactly with the climatology, i.e., dw(xc) = dwc(xc); this corresponds to direct insertion of wc into the model at climatological frequencies. The implementation of spectral nudging in the time domain corresponds to replacing the nudging term by a time convolution of the nudges with a sequence of weights. We denote such a convolved quantity by angle brackets and thus obtain the following time domain form of spectral nudging: wt ¼ wft þ chwct  wft  Uðwct1  wt1 Þi

ð16Þ

Note the main difference between conventional and spectral nudging is that we time-filter the nudges. A critical issue is the design of a computationally efficient and effective filter with zero-phase lag at climatological frequencies. Such a filter is described in Appendix A. The main filter parameter to be specified is the width of the frequency bands over which the nudging occurs. This is controlled by the parameter j (see Appendix A). The reciprocal of j can be interpreted as the time taken, in model time steps, for the filter to spin up. Thus the smaller j, the narrower the width of the nudged bands and the longer the time required for the model to reach a statistical steady state. The other parameter required to implement spectral nudging is c (see (14)). This parameter controls the strength of the nudging. If it is unity, then spectral nudging reduces to direct insertion of the observed climatology at climatological frequencies, as mentioned above. The effect of spectrally nudging the idealized model is illustrated in Fig. 4. For this example jCj = 1 at the climatological frequencies of 0 and 1 cycle per year and j = 103 radians per time step (corresponding to a filter spin-up time of 1000 days and frequency bands of width of 0.2 cycles per year). Note the eddy variability has a peak frequency of 11.4 cycles per year and should not be affected by the spectral nudging. The left panel

116

K.R. Thompson et al. / Ocean Modelling 13 (2006) 109–125

Fig. 4. Spectral nudging of the idealized model. The left panel shows the response of the stream function to the forcing defined in Fig. 1 and with spectral nudging applied according to (16) with c = 1 at 0 and 1 cycle per year, j1 = 1000 d. The middle panel shows the climatology to which the model is nudged (same as Fig. 3). The right panel is the difference between the left and middle panels. Note that it agrees closely with the eddy component of the prognostic run (rightmost panel of Fig. 1) illustrating that spectral nudging does not affect the evolution of the eddy variability. Other model settings are identical to those defined in Fig. 1.

shows the spectrally nudged solution. The middle panel shows the climatology towards which the model is being nudged. The right panel shows the difference between the spectrally-nudged solution and the observed climatology (i.e. the left panel is the sum of the other two panels). Note the difference (the right panel) corresponds very closely to the eddy component in the prognostic run (right panel of Fig. 1). We have therefore shown, for this example at least, that the spectrally nudged solution is essentially the sum of the observed climatology and the eddy component; spectral nudging has had a negligible influence on the evolution of the eddy component as required. The effect of seasonal nudging in the frequency domain is shown in Fig. 2. It is clear that the effect of spectral nudging is confined to two frequency bands centered on 0 and 1 cycle per year; outside of these bands the effect of spectral nudging is effectively zero. The width of the nudged bands is set by j. In this example we have taken a small value of j (corresponding to 0.2 cycles per year) leading to narrow frequency bands over which the model is nudged. It is straightforward to spatially smooth the nudges in addition to applying the temporal filter. The purpose of spatial smoothing is to ensure that the model is nudged only at space scales in the climatology which, in general, will be longer than the dominant scales in the instantaneous model state. Combined with time-filtering, spatial smoothing will help make the eddy variability more transparent to the nudging. Simple spatial filtering corresponds to premultiplying the nudges by an n · n smoothing matrix, S: wt ¼ wft þ cShwct  wft  Uðwct1  wt1 Þi

ð17Þ

K.R. Thompson et al. / Ocean Modelling 13 (2006) 109–125

117

Conceptually this is straightforward and corresponds to nudging the model in specific frequency-wavenumber bands. We have not applied spatial filtering to the idealized model but will apply it to the fully nonlinear, eddy-permitting ocean model discussed in the next section. 5. Spectral nudging of a North Atlantic model The ocean model used to illustrate spectral nudging is Version 2.0 of the Parallel Ocean Program (POP, see Smith et al., 1992). It is a z-coordinate, hydrostatic general ocean circulation model with an implicit free-surface. The model grid covers the North Atlantic basin between 7 N and 67 N. It has 23 vertical levels and a horizontal resolution of 1/3 in longitude and 1/3 cos(/) in latitude (denoted by /). The grid boxes are thus approximately square and range in size from 37 km to 14 km between the southern and northern boundaries. The model is therefore eddy permitting but not eddy resolving over the whole domain. The vertical mixing is parameterized by the KPP scheme of Large et al. (1994) and the horizontal mixing is biharmonic with a tracer mixing coefficient of 5 · 1018 cm4 s1 and a momentum coefficient of 2 · 1019 cm4 s1. The northern and southern boundaries of the model domain are closed. Temperature and salinity are relaxed to the monthly climatology of Yashayaev (see www.mar.dfo-mpo.gc.ca/science/ocean/woce/climatology/naclimatology.htm) within sponge layers with a meridional extension of 3 degrees. The relaxation coefficient is 1/5 d1 at the boundaries of the model domain and smoothly decreases to 0 as the interior is approached. The model was first spun up for 30 years, driven by climatological monthly mean surface fluxes calculated from the COADS data set and spectrally nudged toward the observed seasonal temperature and salinity climatology of Yashayaev. The model was then run several times, with different nudging parameters, for the period 1991–1999 using (i) surface fluxes computed from the model’s sea surface temperature and atmospheric daily mean parameters from the NCEP reanalysis (Kalnay et al., 1996), and (ii) surface water fluxes calculated from the monthly mean precipitation data of Xie and Arkin (1997). For the ‘prognostic run’ the model was not nudged (c = 0), leaving its temperature and salinity fields to evolve freely over the period 1991–1999. For the ‘nudged run’, the model was spectrally nudged toward the observed temperature and salinity climatology using the methodology described in the previous section and Appendix B. Results from the two runs are described and compared below. On examining results from the prognostic run it soon became clear that the model’s mean state was drifting away from the observed climatology and the eddy activity was becoming increasingly unrealistic. This is illustrated in Fig. 5 which shows the standard deviation of the model’s sea level over two 4-year sub-periods. Overall the sea level variances, which are a measure of eddy activity, are too weak and penetrate too far to the east. The mean barotropic stream function (not shown) also indicated that the Gulf Stream’s transport was weakening with time and its path was too far north after leaving the coast at Cape Hatteras. (The poor performance of the model run in prognostic mode was the main motivation for the development of the spectral nudging technique described in this paper.) For the nudged run, both the mean and the annual cycle (xc = 0, 1 cycles per year) were spectrally nudged toward the observed temperature and salinity climatology. The strength of the nudging coefficient c was 1 everywhere except in the top 800 m of the Gulf of Mexico where it was set to zero. (Without this modification the model would not shed Loop Current Eddies. The reason is probably that nudging at a period of 1 year interferes with the eddy shedding process which has a similar time scale.) The filter parameter j1 decreased smoothly through time to the asymptotic value of 5 years. The nudges were spatially smoothed using a fourthorder Butterworth filter. The characteristics of the spatial filter were as follows. Below 50N the filter cutoff was 300 km. North of this latitude the cutoff decreased linearly to 100 km at 60N above which it remained constant at 100 km. The nudges were set to zero over land; after spatial filtering this means that the nudges are attenuated within a few hundred km of the coast. The model temperature at a depth of 87.5 m in the vicinity of the Gulf Stream region is plotted in Fig. 6. Note the model temperature does not drift and its long-term seasonal oscillation follows closely that of the observed climatology (see lower panel). Note also that the filter allows the seasonal cycle to vary from year to year and that the spatial filtering does not require the local seasonal cycle to agree exactly with that of the climatology.

118

K.R. Thompson et al. / Ocean Modelling 13 (2006) 109–125

Fig. 5. Standard deviation of sea level (in cm) from the prognostic run (i.e. no restoring or spectral nudging towards observed temperature and salinity climatologies). The prognostic model integration started in 1991 after a 30 year spin up. The left (right) panel shows the standard deviation over the first (last) part of the model integration. Note the standard deviation generally weakens through time in the vicinity of the Gulf Stream, and also penetrates more to the east.

25

20 °C 15 λ=55W 10 1991

Observed Temperature Modelled Temperature

φ=40N

1992

1993

1994

1995

1996

1997

1998

1999

2000

25

20 °C 15 λ=55W 10 1991

Observed Mean and Annual Cycle Modelled Mean and Annual Cycle

φ=40N

1992

1993

1994

1995

1996

1997

1998

1999

2000

Year

Fig. 6. Time series of observed and modelled temperature at a location in the Northwest Atlantic (40N, 55W and a depth of 87.5 m). The upper panel shows the monthly climatology of Yashayaev (black line) and the modelled water temperature (red line) from the spectrally nudged, North Atlantic model. The observed climatology has a period of 1 year; intramonthly values were estimated by linear interpolation. The modelled temperatures are plotted every 5 d. The lower panel shows the result of fitting the sum of a mean and annual cycle to the observed and modelled temperatures plotted in the upper panel.

The mean sea surface topography, and a snapshot of the model’s sea level, are shown in Fig. 7. The mean clearly shows the Gulf Stream leaving the coast at Cape Hatteras as observed. The subpolar gyre is also evident; the Labrador Current can be seen just offshore of the 1000 m isobath. The snapshot shows a convoluted Gulf Stream path after the Gulf Stream leaves the coast at Cape Hatteras. Well defined cold core rings are evident south of the Stream. Persistent warm core rings are less common in this run; we attribute this to the relatively low horizontal resolution in the Gulf Stream region. The snapshot also shows the Loop Current entering the Gulf of Mexico along with two Loop Current Eddies.

K.R. Thompson et al. / Ocean Modelling 13 (2006) 109–125

119

Fig. 7. The mean sea surface topography, and a snapshot of sea level, from the spectrally nudged ocean model. The mean (left panel) is for the period 1991–1999 inclusive and the snapshot (right panel) is for 15 October, 1999. The black lines show the 1000 m and 4000 m depth contours.

ST. JOHN’S, NFLD.

YARMOUTH

SEAVEY ISLAND

CAPE MAY

CHARLESTON I

DAYTONA BEACH Sea level correction, ηc 20

20

0

0

KEY WEST –20 Jan

Apr

Jul

Oct

–20 Jan

Apr

Jul

Oct

Fig. 8. The seasonal oscillation of sea level observed by 7 coastal tide gauges (black lines) and predicted by the spectrally nudged ocean model (red lines). For each time series plot, the x-axis is time in months and the y-axis is the sea level variation about the long term mean (with the same scale for all plots). The ordering of the subplots follows the latitude of the tide gauges (locations shown by the 7 squares on the map). The observed oscillation was calculated from all monthly sea level data published by the Permanent Sea Level Service for each gauge. (All records were longer than 20 years.) The predicted oscillation is given by gm + gp + gc as explained in the text. The spatially uniform correction, gc in cm, is plotted in the lower left panel.

120

K.R. Thompson et al. / Ocean Modelling 13 (2006) 109–125

To show that the seasonal changes predicted by the model are reasonable we have calculated the observed seasonal oscillation of sea level (go) at 7 representative tide gauges from the western boundary of the North Atlantic (Fig. 8). They were based on monthly means obtained from the Revised Local Reference data directories of the Permanent Service for Sea Level. All monthly mean sea levels for a given month were averaged and then the long term mean was subtracted. One noticeable feature is the large seasonal oscillation at the two gauges in the South Atlantic Bight; this signal has long been attributed to seasonal changes in the Florida Current. Before comparing the model’s sea level (gm) to go, two corrections had to be made. First the inverted barometer effect (gp, see Gill, 1982), which associates a 1 cm increase in sea level with a 1 mb decrease in local air pressure, was added to gm. The second correction accounts for the fact that the model has closed lateral boundaries and conserves volume rather than mass. In such a model, sea level cannot change through mass exchange with the rest of the world’s oceans, or with net heating. In reality both change sea level through a fast barotropic response. The result is that, on seasonal time scales, the associated sea level correction is approximately uniform over the model domain. We made this correction by requiring gm + gp to agree with go at an island station in the North Atlantic where local effects are expected to be small. (We chose St. Georges, Bermuda located at 32.4 N, 64.7 W. About 62 years of monthly data between 1932 and 1999 were available for this station.) The spatially uniform correction, gc = go  gm  gp, had a range of 11 cm (see lower left panel of Fig. 8). The observed and predicted (gm + gp + gc) seasonal oscillations of sea level at the 7 tide gauges are plotted in Fig. 8. Overall the agreement is excellent. Note the good agreement in the South Atlantic Bight which we take as evidence that some of the gross features of the seasonal circulation inferred by the spectrally nudged model are realistic. To further test the realism of the spectrally nudged results we calculated the standard deviation of gridded sea level fields based on delayed-mode altimeter data from the TOPEX/Poseidon, Jason-1, ERS1/2 and ENVISAT satellite missions. The fields covered the period 1993–2001 inclusive and were generated and processed by the Space Oceanography Division of CLS (Collection Localisation Satellites) located in Toulouse, France. (The work was carried out as part of the Environment and Climate EU ENACT project, EVK2-CT200100117, with support from CNES, the French Space agency.) Corrections were made to the along-track data by CLS to suppress residual orbit errors, tides and the inverse barometer effect (Le Traon et al., 1998). The

Fig. 9. Standard deviation (in cm) of sea level measured by altimeters (left panel) and calculated from the spectrally nudged model output (right panel). The altimeter data is for the period 1991–2001. The model standard deviation is based on model output for the period 1991– 1999.

K.R. Thompson et al. / Ocean Modelling 13 (2006) 109–125

121

Fig. 10. Standard deviation (in cm) of sea level measured by altimeters (black line) and calculated from spectrally nudged model output (red line). The meridional sections of standard deviation are taken from the fields plotted in Fig. 9. The spacing between the meridional lines corresponds to 50 cm in standard deviation.

spatial mapping was performed by CLS on the along-track data using an optimal interpolation procedure (Ducet et al., 2000). The standard deviation of the altimeter data is shown in Fig. 9. Maps of this type have been published and discussed extensively in the literature (e.g. Ducet et al., 2000), often with reference to the validation of high resolution ocean models (e.g. Ducet and Le Traon, 2001). We simply note here that the map clearly indicates the region of intense mesoscale variability in the vicinity of the Gulf Stream. The standard deviation of the model sea level (gm) for the period 1991–1999 is shown in the right panel. The semi-quantitative agreement between the two estimates of sea level standard deviation is encouraging and suggests that the mesoscale variability in the model has about the right spatial distribution and amplitude, at least as far as changes in sea level are concerned. A more quantitative comparison of the altimeter and model’s standard deviations is shown by the 9 meridional sections in Fig. 10. This figure shows that the quantitative agreement is reasonable, including the local maxima in the Caribbean and Gulf of Mexico. By way of contrast, the standard deviation of sea level from the prognostic run is too low and spreads too far to the east (compare Figs. 5 and 9). 6. Summary and discussion We have described a computationally efficient and effective scheme for suppressing bias and drift in the temperature and salinity of eddy resolving ocean models. The scheme is based on the simple idea of nudging the ocean model towards observed seasonal climatologies in frequency and wave-number bands where the climatology has energy. Outside of these bands the model evolves prognostically. To test the scheme we compared predictions from a spectrally nudged, eddy permitting model of the North Atlantic against sea levels measured by coastal tide gauges and altimeters. The comparisons were generally encouraging: the model captured the main features in the standard deviation of altimeter observations, and also reproduced the seasonal oscillation of sea level along the Atlantic coast of Canada and the US, including a strong signal in the South Atlantic Bight which is attributed to seasonal changes in the Florida Current. Overall the results suggest that the gross features of the spectrally nudged model, in both nudged and un-nudged bands, are reasonable and sustained over the decadal scale integrations performed here.

122

K.R. Thompson et al. / Ocean Modelling 13 (2006) 109–125

It is important to note that comparison of observations and spectrally nudged output also identified some significant discrepancies. First, the model’s seasonal cycle in sea level in the southern Gulf of Mexico and along the Antilles do not agree well with sea level observations. This suggests the model may not have correctly partitioned the northward flow offshore of the Antilles and through the Caribbean. Another discrepancy was that the standard deviation of model output in the Gulf Stream region was not as peaked, or as continuous along the mean Gulf Stream path, as the altimeter observations. Both discrepancies may be related to the horizontal resolution of the model which is best characterized as eddy permitting rather than eddy resolving. We are presently implementing spectral nudging in higher resolution models of the North Atlantic. Another possible reason for the discrepancies between observed and predicted sea level could be poor lateral boundary conditions. This speculation is based in part on examination of maps of the temperature and salinity nudges at various depths; they suggest large scale distributions that may be traced back to both the northern and southern boundaries. This raises the important point that although spectral nudging can eliminate bias and drift it does so at the expense of adding artificial sources and sinks of heat and salt. Although this may not be important from the perspective of operational modellers interested in forecasts with lead times of days to seasons, it does suggest that the model and its forcing are incorrect in some fundamental ways. We suggest that the nudges, which may be interpreted as flux divergences in the governing equations for the evolution of temperature and salinity, may be useful in diagnosing problems with the model and forcing which, when fixed, will lead to a more physically realistic ocean model with weaker nudging. This is consistent with the point of view expressed by Fox and Haines (2003) in their study of water mass transformations in a global ocean model, and the use of data assimilation in the diagnosis of fundamental problems with the mixing parameterizations, and fluxes through the open boundaries, of the model. Spectral nudging requires the specification of three main parameters. j controls the width of the frequency bands centered on the climatological frequencies (0, 1 cycle per year and possibly its harmonics) to be nudged. Reducing j allows more prognostic variability but also increases the time taken for the time filter to spin up. The second parameter is the strength of the nudging at climatological frequencies. In the present study we have taken it to be unity over most of the domain which, in the linear case, corresponds to direct insertion of the climatological data at the climatological frequencies. The third parameter is the wavenumber cut-off of the spatial smoothing of the nudges. It is not possible to give optimal values for these parameters; the choice will be application dependent. The difficulty in specifying these parameters is akin to the difficulty faced when specifying the covariance structure of observation and model errors. Spectral nudging has potential applications beyond hindcasting and forecasting deep ocean variability. For example it may be useful in shelf modeling where deep ocean eddies would be replaced by diurnal and semidiurnal tides and atmospherically-generated coastal trapped waves. Variations in density stratification are known to affect the dispersion relationship of coastal trapped waves, and also the vertical structure of currents. Thus maintaining the seasonally-varying background density structure by spectral nudging might lead to a more realistic simulation of shelf variability. Another potential application is in coupled atmosphere– ocean modelling where a spectrally-nudged deep ocean model will reduce bias and drift in the ocean component and, through surface fluxes, possibly the atmospheric component of the system as well. Work is underway on both of these applications. Appendix A. A simple, recursive climatological filter with zero phase shift The implementation of spectral nudging requires a computationally efficient filter that passes energy at specific climatological frequencies with zero phase shift. To describe such a filter, let x1, x2, . . . , xt denote a sequence of equispaced scalar values to be filtered. Let st denote an auxiliary state vector that evolves according to st ¼ ðI  KH ÞMst1 þ Kxt where K ¼ j½1 2 0 2 0T H ¼ ½1 1 0 1 0

K.R. Thompson et al. / Ocean Modelling 13 (2006) 109–125

and

2

1

6 60 6 M ¼6 60 6 40 0

0

0

0

0

cosðx1 Þ sinðx1 Þ

 sinðx1 Þ cosðx1 Þ

0 0

0 0

0 0

cosðx2 Þ sinðx2 Þ

123

3

7 7 7 7 7 7  sinðx2 Þ 5 cosðx2 Þ 0 0

where x1 and x2 are the climatological frequencies (henceforth xc) of one and two cycles per year respectively. The filtered value of xt, based on values up to and including time t, is given by hxt i ¼ Hst To describe the frequency response of the filter, assume xt = eixt. The complex amplitude of the filtered response, i.e., the transfer function of the filter, is 1

C ¼ H ½I  ðI  KH ÞMeix  K

ðA:1Þ 1

It is straightforward to show that if j  1 then C approaches [1  i(x  xc)/j] as x approaches xc. Thus at the climatological frequencies, the phase of the filter is zero and the amplitude is unity. It is also straightforward to show the width of the pass band of the filter, p defined as the frequency range over which xt is attenffiffiffi uated by less than one half, is given approximately by 2 3j in radians per model time step. Thus the smaller j, the narrower the pass band of the filter. The spin-up of the filter is also controlled by j. It is straightforward to show that as j ! 0+ the eigenvalues of (I  KH)M tend to 1  j in absolute value. This means that the spin-up time of the filter is of order j1 in model timesteps. For linear ocean models, such as the one described in Section 2, the transfer function of the spectrally nudged model is given by the following frequency dependent generalization of (11) for x5xc:

1

Gain

0.8 0.6 0.4 0.2

Phase, in degrees

0

0

0.5

1

1.5

2

2.5

3

0

0.5

1

1.5 Frequency, cycles per year

2

2.5

3

50

0

–50

Fig. 11. Transfer function of the filter used in the spectrally nudged ocean model. The top panel shows the gain (j1  Cj) and the bottom panel shows the phase (arg(1  C), in degrees) where C is the transfer function of the climatological filter defined by (A.1). The shaded areas mark the regions that are within j of the nudged frequencies (0 and 1 cycle per year in this example). We took j1 = 1000 d.

124

K.R. Thompson et al. / Ocean Modelling 13 (2006) 109–125

T s ¼ ð1  cCÞT u

ðA:2Þ

where Tu is the transfer function of the un-nudged model. A plot of the gain and phase of 1  C for j1 = 1000 d is shown in Fig. 11. Note that outside of the nudged bands (i.e. the shaded areas) the gain and phase of 1  C are approximately 1 and 0 respectively showing that spectral nudging has little effect there (i.e. Tn  Tu). The above equations provide a simple, recursive way of filtering any input time series. In the present context the input is a function of the difference between the observed climatology and the model prediction. The filter passes variability in the vicinity of the specified climatological frequencies; the width of the pass band is controlled by j. The strength of the nudging is controlled by the parameter c. If c = 1 the longterm climatology of the model will agree with the observed climatology at x = xc; if c = 0, no nudging is applied and the model is prognostic. Several extensions are possible. For example it is straightforward to add additional harmonics to the annual cycle and allow both j and c to vary with time, location and xc (e.g. nudge the mean more strongly than the annual cycle). Appendix B. Implementation of spectral nudging in a general circulation model This appendix outlines the implementation of spectral nudging in the Parallel Ocean Program (POP). We begin with the analogue of (17) written for a single tracer x that represents either temperature or salinity in the POP model: xtþ1 ¼ xftþ1 þ cShxctþ1  xftþ1  /ðxct1 ; xct ; ut Þ þ /ðxt1 ; xt ; ut Þi

ðB:1Þ

c

The terminology is the same as (17) except that (i) x represents the observed tracer climatology at the grid point and time of interest, and (ii) the arguments of the dynamics operator / include tracer values from the previous two time levels, and the velocity from the central time level, ut, as required for leapfrog time-stepping of the advection–diffusion equation. The last term on the right side of (B.1) is the spectral nudging term; the angled brackets represent the recursive filter described in Appendix A and S represents the spatial smoothing operator. In the POP model run discussed in Section 5, we used the climatological filter described in Appendix A at each time step, and a low-pass Butterworth filter in both zonal and meridional directions at each model level. We next simplify (B.1) by replacing /ðxct1 ; xct ; ut Þ by axct , and /(xt1, xt, ut) by axt, where a is a scalar to be specified. The unforced advection–diffusion equation for the tracers is thus approximated by a simple scalar multiplication within the nudging term, obviously a crude approximation. This gives xtþ1 ¼ xftþ1 þ cShxctþ1  xftþ1  aðxct  xt Þi

ðB:2Þ

Note from (B.2) that setting a = 0 leads to a frequency-wavenumber generalization of conventional nudging as defined by (7). A non-zero value of a provides a linearization of the nudging approach corresponding to (17). In the application described in Section 5 we took a = 0.99 (for an equivalent time step of 1 day). The steps required to implement spectral nudging in a complex ocean model are thus as follows: (1) Update all prognostic variables to time t + 1 using the ocean model. This gives the one-step-ahead forecast, xftþ1 . (2) Calculate xctþ1  xftþ1  aðxct  xt Þ for each T and S grid point and filter in time to obtain hxctþ1  xftþ1  aðxct  xt Þi. (3) Spatially smooth the filtered nudges to obtain Shxctþ1  xftþ1  aðxct  xt Þi. (4) Update the one-step-ahead forecast, i.e., calculate xtþ1 ¼ xftþ1 þ cShxctþ1  xftþ1  aðxct  xt Þi. (5) Save xctþ1  xtþ1 for step 2 on the next time level. The above steps must be carried out for each tracer grid point. For the present North Atlantic application, which is based on 1/3 degree horizontal resolution and 23 levels, the computational cost of spectral nudging is about 20% of the total cost. The memory requirements are only a few percent of the total memory requirements of the model; they come from the vector of length 5 (i.e. st in Appendix A) that must be carried for each grid variable to be nudged.

K.R. Thompson et al. / Ocean Modelling 13 (2006) 109–125

125

References Bleck, R., 2002. An oceanic general circulation model framed in hybrid isopycnic-Cartesian coordinates. Ocean Model. 37, 55–88. Cooper, M., Haines, K., 1996. Altimetric assimilation with water property conservation. J. Geophys. Res. 101, 1059–1077. Daley, R., 1994. Atmospheric Data Analysis. Cambridge University Press. Ducet, N., Le Traon, P.Y., 2001. A comparison of surface eddy kinetic energy and Reynolds stresses in the Gulf Stream and the Kuroshio current systems from merged TOPEX/Poseidon and ERS-1/2 altimetric data. J. Geophys. Res. 106 (C8), 16,603–16,662. Ducet, N., Le Traon, P.Y., Reverdin, G., 2000. Global high-resolution mapping of ocean circulation from TOPEX/Poseidon and ERS-1 and 2. J. Geophys. Res. 105 (C8), 19,477–19,498. Evensen, G., 1994. Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics. J. Geophys. Res. 99 (5), 10143–10162. Fox, A.D., Haines, K., 2003. Interpretation of water mass transformations diagnosed from data assimilation. J. Phys. Oceanogr. 33 (3), 485–498. Gill, A.E., 1982. Atmosphere–Ocean Dynamics. Academic Press, New York. Halliwell, G.R., 2004. Evaluation of vertical coordinate and vertical mixing algorithms in the HYbrid-Coordinate Ocean Model HYCOM. Ocean Model. 7, 285–322. Huddleston, M.R., Bell, M.J., Martin, M.J., Nichols, N.K., 2004. Assessment of wind-stress errors using bias corrected ocean data assimilation. Quart. J. Roy. Meteorol. Soc. 130, 853–871. Jazwinski, A.H., 1970. Stochastic Processes and Filtering Theory. Academic Press, New York. Kalnay, E., Kanamitsu, M., Kistler, R., Collins, W., Deaven, D., Gandin, L., Iredell, M., Saha, S., White, G., Woollen, J., Zhu, Y., Chelliah, M., Ebisuzaki, W., Higgins, W., Janowiak, J., Mo, K.C., Ropelewski, C., Wang, J., Leetmaa, A., Reynolds, R., Jenne, R., Joseph, D., 1996. The NCEP/NCAR 40-year reanalysis project. Bull. Am. Meteorol. Soc. 77, 437–471. Kitagawa, G., 1998. A self-organizing state-space model. J. Am. Statist. Assoc. 93 (443), 1203–1215. Large, W.G., McWilliams, J.C., Doney, S.C., 1994. Oceanic vertical mixing: a review and a model with a non-local K-profile boundary layer parameterization. Rev. Geophys. 32, 363–403. Le Traon, P.-Y., Nadal, F., Ducet, N., 1998. An improved mapping method of multisatellite altimeter data. J. Atmos. Ocean. Technol. 15, 522–534. Longuet-Higgins, M.S., 1964. Planetary waves on a rotating sphere. Proc. Roy. Soc. London A 279, 446–473. Lorenc, A., 1986. Analysis methods for numerical weather prediction. Quart. J. Roy. Meteorol. Soc. 112, 1177–1194. Pedlosky, J., 1965. A note on the western intensification of the oceanic circulation. J. Mar. Res. 23, 207–209. Pham, D.T., Verron, J., Roubaud, M.C., 1998. Singular evolutive extended Kalman filter with EOF initialization for data assimilation in oceanography. J. Mar. Syst. 16 (3–4), 323–340. Priestley, M.B., 1981. Spectral Analysis and Time Series, vol. 1. Academic Press. Sarkisyan, A.S., 1977. The diagnostic calculations of a large-scale oceanic circulation. In: Goldberg, E.D., Mccave, I.N., O’Brien, J.J. (Eds.), The Sea, vol. 6. Wiley-Interscience Publication, New York, pp. 363–458. Sarmiento, J.L., Bryan, K., 1982. An ocean transport model for the North Atlantic. J. Geophys. Res. 87, 394–408. Smith, R.D., Dukowicz, J.K., Malone, R.C., 1992. Parallel ocean general circulation modeling. Physica D 60, 38–61. Thacker, W.C., Esenkov, O.E., 2002. Assimilating XBT Data into HY-COM. J. Atmos. Ocean. Technol. 19 (5), 709–724. von Storch, H., Langenberg, H., Feser, F., 2000. A spectral nudging technique for dynamical downscaling purposes. Monthly Weather Rev. 128, 3664–3673. Woodgate, R.A., Killworth, P.D., 1997. The effects of assimilation on the physics of an ocean model. Part 1: Theoretical model and barotropic results. J. Atmos. Ocean Technol. 14, 897–909. Xie, P., Arkin, P.A., 1997. Global precipitation: a 17-year monthly analysis based on gauge observations, satellite estimates, and numerical model outputs. Bull. Am. Meteorol. Soc. 78, 2539–2558.