A reduced-dynamics variational approach for the assimilation of altimeter data into eddy-resolving ocean models

A reduced-dynamics variational approach for the assimilation of altimeter data into eddy-resolving ocean models

Ocean Modelling 27 (2009) 215–229 Contents lists available at ScienceDirect Ocean Modelling journal homepage: www.elsevier.com/locate/ocemod A redu...

3MB Sizes 0 Downloads 57 Views

Ocean Modelling 27 (2009) 215–229

Contents lists available at ScienceDirect

Ocean Modelling journal homepage: www.elsevier.com/locate/ocemod

A reduced-dynamics variational approach for the assimilation of altimeter data into eddy-resolving ocean models Peng Yu 1, Steven L. Morey *, James J. O’Brien Center for Ocean – Atmospheric Prediction Studies, The Florida State University, Tallahassee, FL 32306-2840, USA

a r t i c l e

i n f o

Article history: Received 4 January 2008 Received in revised form 11 January 2009 Accepted 16 January 2009 Available online xxxx Keywords: Ocean modeling Data assimilation Variational adjoint methods

a b s t r a c t A new method of assimilating sea surface height (SSH) data into ocean models is introduced and tested. Many features observable by satellite altimetry are approximated by the first baroclinic mode over much of the ocean, especially in the lower (but non-equatorial) and mid latitude regions. Based on this dynamical trait, a reduced-dynamics adjoint technique is developed and implemented with a three-dimensional model using vertical normal mode decomposition. To reduce the complexity of the variational data assimilation problem, the adjoint equations are based on a one-active-layer reduced-gravity model, which approximates the first baroclinic mode, as opposed to the full three-dimensional model equations. The reduced dimensionality of the adjoint model leads to lower computational cost than a traditional variational data assimilation algorithm. The technique is applicable to regions of the ocean where the SSH variability is dominated by the first baroclinic mode. The adjustment of the first baroclinic mode model fields dynamically transfers the SSH information to the deep ocean layers. The technique is developed in a modular fashion that can be readily implemented with many three-dimensional ocean models. For this study, the method is tested with the Navy Coastal Ocean Model (NCOM) configured to simulate the Gulf of Mexico. Ó 2009 Elsevier Ltd. All rights reserved.

1. Introduction Measurements of the ocean environment are difficult and expensive, especially within the deep ocean. As a result, oceanic data are usually sparse and non-uniformly sampled in time and space. The launch of many ocean observing satellites over recent decades has yielded a vast amount of data for the ocean’s surface. In order to deepen and broaden our understanding of ocean circulation, it is very important to optimize the use of these expanded valuable datasets. Data assimilation is a powerful tool for extracting the maximum amount of information from observations, blending present observations with the theoretical knowledge from past observations (Le Dimet and Navon, 1988; Gill and Malanotte-Rizzoli, 1991; Wunsch 1996), and has been extensively used in numerical oceanic modeling. In order to assimilate observations into numerical ocean models, there exist a variety of different methods, among which the variational method (Le Dimet and Talagrand, 1986; Courtier and Talagrand, 1990; Yu and O’Brien, 1991; Chua and Bennett, 2001; Bennett, 2002; Wunsch and Heimbach, 2007) and Kalman filter (Ghil et al., 1981; Evensen, 1994; Fukumori and Malanotte-Rizzoli, 1995; Fukumori, 2002) are

* Corresponding author. Tel.: +1 850 644 0345; fax: +1 850 644 4841. E-mail address: [email protected] (S.L. Morey). 1 Present address: Jet Propulsion Laboratory, California Institute of Technology, M/S 300-323, 4800 Oak Grove Drive, Pasadena, CA 91109, USA. 1463-5003/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.ocemod.2009.01.006

the two most advanced approaches. Variational methods estimate a system state that fits a dynamical model and available observations in a least-square sense. The four-dimensional variational data assimilation (4D-VAR) minimizes a cost function, the distance between the observations and their model counterpart, by adjusting some chosen model parameters (control variables), such as initial conditions, while the set of model dynamics acts as a strong constraint. It is very resource-demanding to solve a complicated 4D-VAR problem in meteorology and oceanography, especially with a high-resolution eddy-permitting model. Reduced-space approximations, resulting in fewer degrees of freedom than those formally in the models, have been applied to reduce the computational cost. Primarily, these techniques reduce the order of the problem by projecting the state space onto a linear subspace. One commonly used approach is to construct a low dimension subspace based on the first few empirical orthogonal functions (EOF) (Blayo et al., 2003; Robert et al., 2005; Hoteit, 2008). Hoteit et al. (2005) introduced a decomposition of the velocities at the open boundaries to deal with very strong sensitivities of the model sea surface height (SSH) to the barotropic component of the inflow. Particularly, Kohl and Willebrand (2002) improved the state estimate of a dynamical system by the assimilation of statistical moments of the model state. The purpose of this study is to develop a computationally efficient variational assimilation approach by using an adjoint model

216

P. Yu et al. / Ocean Modelling 27 (2009) 215–229

with a reduced state space. Different from the above-mentioned techniques, the subspace in this new method is constructed based on the local ocean dynamics. The adjoint model is simplified by approximating the ocean model variable fields by their first baroclinic mode counterparts. Vertical normal mode decomposition and reconstruction is employed to link the different forward and backward models with different dynamics. Because of the near-surface intensification of baroclinic modes, sea level perturbations measured by altimeters primarily reflect the first baroclinic mode, and thus the motion of the main thermocline (Wunsch, 1996). The application of the method here assimilates satellite altimeter data and provides a mechanism for transferring the information vertically from the ocean surface to deep layers, leading to an improved estimate of the ocean state throughout the water column. Complementing the assimilation method is an application of an interpolation algorithm to map the inhomogeneously sampled SSH to horizontally propagate the information from observation locations to data void areas and thus expedite the assimilation procedure. Several assumptions are made in this reduced-dynamics variational data assimilation technique, which lead to a more efficient way of assimilating the observations into the model system at a cost of an approximated gradient, and hence a compromised solution. Three numerical experiments are conducted with a Gulf of Mexico (GoM) configuration of the Navy Coastal Ocean Model (NCOM) for testing purposes. The first experiment assimilates the full SSH field from a numerical model solution as ‘‘observations” into a model with a perturbed initial condition in attempt to recover the original initial conditions. This experiment is designed to test the reduced-dynamics assimilation method, show that it can handle a poorly chosen initial guess of the ocean state, and demonstrate the assimilation of surface fields into the ocean model’s three-dimensional state variables. A similar experiment, but with the inclusion of a mapping technique, is tested in the second trial. This experiment simulates a more realistic scenario by sampling the unperturbed model SSH as observations in a similar manner as the satellite altimeters (along satellite ground tracks). The sampled SSH data are mapped to a regular grid offline using a mapping technique before the assimilation procedures. The mapping technique horizontally brings the information from observation locations to the native model grid, and thus expedites the horizontal propagation of the observation information in the assimilation. The last experiment is a real-world implementation of the technique using observations from TOPEX/Poseidon (T/P) and Jason-1 satellite altimeters.

2. Methodology and data The assimilation system consists of several components, combined in a modular fashion that can be applied to a number of three-dimensional ocean models. For this research, the NCOM is used as the ‘‘forward” model for testing the method. The adjoint, often called the ‘‘backward model” because it is integrated backward in time, is developed from a simplified model, a one-activelayer reduced-gravity model. The forward and backward models are connected through vertical normal mode decomposition and reconstruction. These components form the data assimilation system when iterated within an optimization algorithm (Fig. 1). This new technique is designed to efficiently assimilate satellite altimeter data into the three-dimensional ocean model. The general form of the cost function in 4D-VAR problems is defined as the sum of model background error and the observation error. A well-defined background error term can greatly improve the performance of the system. However, because of the lack of knowledge of the model error, it is usually difficult to define the error covariance matrix. Navon et al. (2005) investigated the impact of the inclusion of the background error on incomplete

observations (observations are more sparse than the model stencil) for 4D-VAR data assimilation. The results suggested that when observations are complete in space, the inclusion of the background error is not important, and may lead to negative impacts in some cases. In this study, a statistical interpolator based on complex EOF analysis is used to interpolate the satellite altimeter data to the model native grid. As a result, the background error term is excluded from the definition of the cost function in this research. 2.1. The reduced-dynamics variational data assimilation method Oceanic state variables can be dynamically decomposed into different vertical normal modes, where each individual mode or the combination of several modes accounts for different phenomena. Previous work (Dewar and Morris, 2000; Dewar and Huang, 2001; Shu and Clarke, 2002) shows that many oceanic features of interest can be approximated by the first few baroclinic modes. It has been suggested (Liu, 1999) that the SSH variability is mainly associated with the first baroclinic mode. Wunsch (1997) inferred that surface kinetic energies are dominated by the first baroclinic mode, and thus to a first order approximation, the altimeter over much of the World Ocean directly reflects the movement of the main thermocline. Siegel et al. (1999) showed that in the subtropical North Atlantic, the first baroclinic mode explained 88.5% of the SSH anomaly variance using hydrographic and satellite altimeter data. A similar analysis conducted within the GoM model domain used for this present research reveals that the first baroclinic mode accounts for around 80% of the total variance. As a result, it is reasonable to reduce the dimensions of the assimilation problem by using a reduced-gravity adjoint model, which closely approximates the first baroclinic mode, leading to an efficient assimilation scheme. This assumption restricts the application of the technique only to SSH observations; however, these data are most commonly used for assimilation into operational ocean models as they have been collected continuously for a long period of time and can yield considerable information regarding the ocean dynamics. The variational assimilation procedure works by minimizing a cost function, the data misfit between the SSH observations and their counterpart model variables. For this research, the assimilation method is applied to a three-dimensional time-dependent forward ocean model that is integrated within the assimilation window. Instead of constructing and running an adjoint model based on the fully three-dimensional forward model as in the conventional adjoint data assimilation technique, an adjoint model based on a one-active-layer reduced-gravity model, which approximates the first baroclinic mode, is developed. The different forward and backward models are dynamically connected through vertical normal mode decomposition, with which the first baroclinic mode information is retrieved from the prognostic variables of the forward model and then used when integrating the adjoint model. After the integration of the approximate adjoint model, the gradient of the cost function with respect to the control variables (velocity and SSH associated with the first baroclinic mode) is obtained as an input to an optimization routine. For this study, the limited memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) method (Liu and Norcedal, 1989), which is of the class of Limited Memory Quasi-Newton (LMQN) methods, is used. The optimization algorithm adjusts the first baroclinic mode velocity and SSH fields of the initial conditions to reduce the magnitude of the cost function. The new values of these fields are then used to update the full forward model initial variables by reconstructing the model fields from the vertical modes (with the newly updated first mode) as described in Section 2.2. This procedure repeats until the decrease of the cost function value of the current iteration over that of the previous one normalized by the current cost function value is less than a prescribed value.

217

P. Yu et al. / Ocean Modelling 27 (2009) 215–229

Ocean Model (0-T)

Output Files

Altimeter Data

Mapping Method

Gridded Data

Vertical Normal Mode Decomposition

ˆ ˆ h) (u,v

Reduced Gravity Adjoint Model (T-0)

REDECED SPACE ADJOINT DATA ASSIMILATION

No

MINIMIZATION PROCEDURE: ITERATE UNTIL CONVERGENCE

Input Files (Init. Conditions)

Vertical Normal Mode Reconstuction

ˆ ˆ h)t=0 (u,v

Yes Forecast (T-T 2) Fig. 1. Flowchart of the reduced-dynamics variational data assimilation technique. The models are run for the assimilation window t = [0, T], and the forward model solution at time t = T can be used to initialize a forecast to some time t2.

2.2. Vertical normal mode decomposition Vertical normal mode decomposition is applied to retrieve the information from the forward ocean model. The vertical normal modes are determined by solving the following partial differential equation of Sturm–Liouville type:

Szz þ

N2 S¼0 gh

ð2:1Þ

with boundary conditions:

S ¼ 0; at z ¼ zbottom ; at z ¼ 0; S  hSz ¼ 0;

ð2:2Þ ð2:3Þ

(see Appendix A for the detailed derivation of (2.1)–(2.3)), where S is a function of z, N is the Brunt-Väisälä frequency, g is the acceleration of gravity, and h is the equivalent depth. Discretization of (2.1)–(2.3) yields an eigenvalue problem:

AS ¼ kS; where A is an n  n tri-diagonal matrix, k is the eigenvalue, and S is the eigenvector. The first eigenmode, called mode zero, represents the barotropic mode and all the other modes are baroclinic. Based on the magnitude of the eigenvalues, the baroclinic modes are ordered as mode one, mode two, etc. The mode number m is the number of zero-crossing(s) of the eigenvector Rm (Fig. 2). The eigenvectors Rm and Sm give the vertical structure of each mode, and the equivalent depth associated with each mode is obtained from the corresponding eigenvalue. Once the vertical normal modes are known, the complete three-dimensional timedependent fields can be reconstructed as M X



U m ðx; y; tÞRm ðzÞ;

ð2:4Þ

V m ðx; y; tÞRm ðzÞ;

ð2:5Þ

m¼0



M X m¼0

p g q0 w¼

¼

M X

gm ðx; y; tÞRm ðzÞ;

ð2:6Þ

m¼0 M X m¼0

W m ðx; y; tÞSm ðzÞ:

ð2:7Þ

Suppose at the end of an iteration of the optimization algorithm, the increments of the first baroclinic mode variables are DU 1 , DV 1 , and Dg1 . Since only the first baroclinic mode part of the field is updated and the separation of variables is linear, the two horizontal velocity components and the pressure field are updated in the form of

uiþ1 ðx; y; z; 0Þ ¼ ui ðx; y; z; 0Þ þ DU 1 ðx; y; 0ÞR1 ðzÞ; iþ1

i

ð2:8Þ

ðx; y; z; 0Þ ¼ v ðx; y; z; 0Þ þ DV 1 ðx; y; 0ÞR1 ðzÞ;

ð2:9Þ

piþ1 ðx; y; z; 0Þ ¼ pi ðx; y; z; 0Þ þ g q0 Dg1 ðx; y; 0ÞR1 ðzÞ;

ð2:10Þ

v

based on (2.4)–(2.6), where i is the iteration number. The update of the density field, and hence temperature and salinity fields, at the model initial time is crucial for the forward model solution, since the density field has the most immediate impact on the ocean surface elevation and circulation patters. Cooper and Haines (1996) described a way to project the update of the SSH field to that of the density field in the assimilation procedures under the assumption that the ocean bottom remains a level of no motion (i.e., the change in weight of the entire water column should compensate for the change in surface pressure observed by the altimeter). Similarly, through the calculation of the vertical profile of the first baroclinic mode, the SSH update can be easily projected to the density profile. Differentiating (2.10) with respect to z, one gets i piþ1 z ðx; y; z; 0Þ ¼ pz ðx; y; z; 0Þ þ g q0 Dg1 ðx; y; 0Þ

dR1 ðzÞ : dz

ð2:11Þ

From (A17), the hydrostatic balance, the density update is

qiþ1 ðx; y; z; 0Þ ¼ qi ðx; y; z; 0Þ  q0 Dg1 ðx; y; 0Þ

dR1 ðzÞ : dz

ð2:12Þ

The updated density needs to be converted to temperature and salinity variables for the forward model. For the application of the method used in this research, a mapping from density onto temperature and salinity is derived from the local temperature and salinity profile of the forward model as follows. At each grid point, the local vertical density, temperature, and salinity profiles averaged over time are obtained from the previous model run and the relation of these three variables is derived using the regression method. The update of the temperature and salinity variables is based on this relation and the density values. A quality control pro-

218

P. Yu et al. / Ocean Modelling 27 (2009) 215–229

Density Profile

First Three Baroclinic Modes 0

200

200

400

400 Depth (m)

Depth (m)

0

600

600

800

800

1000 1024

1025 1026 1027 Density (kg/m3)

1028

1000 −2

0

2 Amplitude

4

6 x 10−3

Fig. 2. The average density profile of the Gulf of Mexico model domain to which vertical normal mode decomposition is applied. Right: Eigenvectors R associated with the first three baroclinic modes. Solid line, mode one; dashed line, mode two; dotted line, mode three.

cedure is conducted after this, since non-realistic values may rarely occur at some locations because of overshooting. An extreme value is replaced with the weighted average of the values from its neighboring points. 2.3. Numerical optimization The efficient implementation of the variational assimilation method depends upon the fast convergence of a large-scale unconstrained minimization algorithm. A review of experiences and details concerning various algorithms can be found in Zou et al. (1993). Since problems in oceanography contain many degrees of freedom, the Limited Memory Quasi-Newton method (LMQN) (Navon and Legler, 1987; Navon et al., 1992) is chosen because of its ability to deal with large-scale problems efficiently. In the LMQN algorithm, the optimal solution is approximated asymptotically by adjusting the value of the control variables iteratively along its gradient or a direction obtained from a combination of the gradients from previous iterations. Various LMQN methods exist, and the L-BFGS method (Liu and Norcedal, 1989) is used in this study.

considers not only the observations within the assimilation cycle but also those prior to the assimilation window. Though this intermediate procedure is not required, the inclusion of it allows a shorter assimilation period, since the technique maps the original observations from satellite ground track observational locations onto a finer spatial and temporal grid and thus reduces the time for the spreading of information to unobserved areas. 2.5. The ocean model The ocean model used to test the reduced-dynamics variational data assimilation technique is the NCOM developed at the Naval Research Laboratory (Martin, 2000; Morey et al., 2003). It uses a hybrid sigma-z level coordinate. The free surface is a sigma level, and the levels below the free surface are a specified number of sigma and z-levels. The NCOM is a fully three-dimensional primitive equation ocean model with the hydrostatic, incompressible, and Boussinesq approximations. The model supports a number of higher order numerical methods. For these simulations, a quasi-thirdorder advection scheme is used with a leapfrog time differencing scheme. Details of the simulations are described in Section 3.1.

2.4. TOPEX/Poseidon and Jason-1 altimeter data Data from the T/P and Jason-1 satellite altimeters are assimilated for purposes of testing and demonstration. The T/P satellite altimeter provided SSH measurements with a root mean square error (RMSE) of less than five centimeters (around 4.2 cm) and the Jason-1 around 3.9 cm. During the period when both satellites were in operation, the spatial coverage of the observations was doubled. Despite the increased spatial coverage provided by the two orbiting altimeters (Fig. 3), there are still large gaps between the ground tracks comparable in size to mesoscale oceanic features. A mapping technique developed by Yu et al. (2004) based on complex empirical orthogonal function analysis (Barnett, 1983; Shriver et al., 1991) is used to produce an SSH data set with a regular temporal and spatial grid for assimilation. The mapping method

3. Experiments and results The reduced-dynamics variational data assimilation method, described in Chapter 2, is applied to an NCOM simulation of the GoM for testing and demonstration. Three numerical experiments are run. The first is an idealized experiment in that the SSH data to be assimilated are from a model run. The data are then assimilated into a model that has started from a different initial state and the goal is to adjust the model initial state toward the original one from which the assimilated SSH data are obtained. The difference between the initial conditions of the control run and the first guess are large, so as to test how the system responds to a bad first guess. In particular, the sizable density perturbations are a good test of

219

P. Yu et al. / Ocean Modelling 27 (2009) 215–229

T/P and Jason 1 ground tracks in the Gulf of Mexico 0m Louisiana

30˚N

MAFLA MAFLAShelf Shelf Texas

LATEX LATEXShelf Shelf

Mexico Mexico East M e xico Sh 20 elf 00

Latitude

West West Florida Florida Shelf Shelf

1000 2000

1000

24˚N

Florida

1000

2000 Campeche CampecheBank Bank Bay of Campeche

Cuba Cuba

Yucatan Peninsula Peninsula

3000

18˚N

4000 98˚W

90˚W Longitude

82˚W

Fig. 3. TOPEX/Poseidon and Jason-1 ground tracks over the Gulf of Mexico model domain (98.15°W–80.60°W, 15.55°N–31.50°N) and topography for the NCOM simulations used in this research. Dashed lines: T/P ground tracks; Solid lines: Jason-1 ground tracks.

the method for projecting the density perturbations onto temperature and salinity profiles. The second experiment is similar to the first, except that the SSH data are sampled along a pattern similar to the T/P and Jason-1 satellite tracks, and the EOF-based mapping technique is applied to construct the gridded SSH fields to be assimilated. In these two idealized experiments, we assume that observation errors are uncorrelated and the error variances are homogeneous. The last experiment is an application of the technique with real data, assimilating T/P and Jason-1 along-track data from January 9, 2004 to January 18, 2004. 3.1. Idealized experiment one The data assimilation technique is tested using the NCOM configured for a domain covering the entire GoM and the northwestern Caribbean (98.15°W–80.60°W, 15.55°N–31.50°N) (Fig. 3). The horizontal resolution is 0.05° in both latitude and longitude and there are 352 grid points zonally and 320 grid points meridionally. It has 20 evenly spaced sigma levels above 100 m depth and 40 zlevels below 100 m with a maximum depth of 4000 m. Open boundaries are found along the eastern edge of the domain in the Caribbean and the Florida Straits. This model is integrated for a period of 11 years from rest forced by wind stress, latent, sensible, and radiative heat fluxes derived from the DaSilva et al. (1994) 0.5°  0.5° analyzed COADS monthly climatology fields. The 1994 World Ocean Atlas (National Oceanic and Atmospheric Administration, 1994) (WOA94) is used to derive the model initial temperature and salinity fields. More details of this model configuration can be found in Morey et al. (2005). The model fields at the end of the 11th year of model integration are used as the initial condition for the first iteration of the assimilation technique (a ‘‘first guess”). The NCOM is then integrated with a different initial condition as a control or truth run. The SSH field at the native model resolution is sampled at

one-day time intervals as the ‘‘observation” fields (Fig. 5a). This experiment assimilates the SSH data from the truth run into the model to adjust it toward the same state as the truth run. This experiment assimilates an SSH observation field into a model that has not been constrained to match the truth fields in any way. Thus, the first guess is quite different from the observation. The purpose is to illustrate how the technique works and to verify its utility given a poor first guess of the initial condition. During each iteration the model is run from an initial condition for a ten-day period and the cost function, defined as in Appendix B, is calculated. Vertical normal mode decomposition is applied to retrieve the first baroclinic mode SSH and horizontal velocity fields as the model control variables and to obtain the eigenfunction corresponding to the first baroclinic mode. The adjoint model is integrated backward in time for 10 days to obtain the gradient of the cost function with respect to the model control variables. The LBFGS method is applied to minimize the cost function and update the first baroclinic mode SSH and horizontal velocity fields at the model initial time. These updated first baroclinic mode variables are used to calculate the update for the full model velocity and pressure fields through vertical normal mode reconstruction. The density, temperature, and salinity are updated based on the algorithm described in Section 2. The whole procedure repeats until a prescribed convergence criterion is met. The variance of the residual SSH field (the difference between the observation and the model counterpart) (Fig. 4a) does not decrease dramatically during the first three iterations, and then drops rapidly in the following several iterations when the optimization algorithm gathers more information from the previous iterations. By iteration 7, it drops to 20% of the value from the first iteration with negligible change following subsequent iterations. The cost function behaves in a similar way as the variance of the residual field. The gradients of the cost function with respect to the control variables are approximated with the first baroclinic mode. The

220

P. Yu et al. / Ocean Modelling 27 (2009) 215–229

at the sea surface to deep layers. Since the adjoint model has only one vertical layer, it takes much less computational time (around 1/30 of that required by the NCOM or its adjoint) than the forward model run resulting in an efficient implementation of a variational data assimilation technique.

(a) Experiment 1 0.05

m2

0.04 0.03

3.2. Idealized experiment two

0.02 0.01 0 1

2

3

4

5 6 Iteration

7

5

6 7 8 Iteration

9

8

9

10

(b) Experiment 2 0.01

m2

0.008 0.006 0.004 0.002 0 1

2

3

4

10 11 12 13

(c) Experiment 3 0.05

m2

0.04 0.03 0.02 0.01 0 1

2

3

4

5

6

7 8 9 10 11 12 13 14 Iteration

Fig. 4. Variance of the residual SSH field (the difference between the observation and the model counterpart) with respect to iteration number from idealized experiment one (a), experiment two (b), and experiment three (c).

inconsistency between the forward and backward model may be the cause of the non-monotonically decreasing cost function. The assimilation technique improves the consistency between the model SSH and the observations and most of the mesoscale variability is captured correctly (Fig. 5). The biggest difference between the first guess and the observational field is the large anticyclonic eddy in the northwestern GoM. It disappears after the assimilation procedure and the SSH field in the area around Cuba improves dramatically. The root mean square error (RMSE, which qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

ðSSHðx; y; tÞ  SSHobs ðx; y; tÞÞ2 where the is defined in the form overbar denotes a time average over the assimilation window) fields (Fig. 5c) supply a better comparison. In order to assess the model/truth fit of other variables, we compared the 18 °C isotherm depth in the model domain. Fig. 6 shows that the main features are much better represented after the assimilation, which indicates that the reconstruction of the fields following assimilation of SSH data into the first baroclinic mode is capable of dealing with rather large density perturbations and transfer the information observed

This experiment is designed to simulate a more realistic scenario in which the SSH is sparsely sampled by satellite altimeters. The model configuration is the same as in experiment one. After the 11-year spin-up period the NCOM is run for an additional 100 days. Rather than assimilating the SSH field sampled at the native model grid, the model SSH data are sub-sampled along locations corresponding to T/P and Jason-1 observation locations every ten days. The mapping technique is applied to interpolate the sub-sampled data (representing satellite altimeter along-track observations) to the native model grid in space and to a one-day interval over time between model year day 90 and day 99. The comparison between the mapped SSH data and the model ‘‘truth” SSH fields (Fig. 7) shows that the gridded SSH constructed from the mapping method represents the original field well. The interpolator introduces errors with standard deviation over the domain of 1 cm. The NCOM model state at the end of the 11-year spin-up period (before the additional 100 days of integration) is taken as the first guess of the assimilation and the reduced dynamics variational data assimilation procedure is applied. Thus, the first guess is a perturbation from the truth as they are separated in time by 90 days. By taking a first guess of the initial condition closer to the truth field than in the previous experiment, we assume that we have a better knowledge of the ocean state at the model initial time. The variance of the residual SSH field reaches 43% of its initial value after iteration six with a rapid reduction of the cost function during the first four iterations (Fig. 4b). The comparison between the truth and the NCOM outputs (Fig. 8a) shows the improvement of the SSH field during the assimilation, e.g., the shape of the large anticyclonic eddy in the northwestern Gulf is closer to the truth, the positions and shapes of the eddies at Bay of Campeche are improved, and the shedding of an anticyclonic eddy at the end of the assimilation window is recovered. The surface velocity field is also dynamically consistent with the SSH field. The SSH RMSE (Fig. 8b) gives a clearer comparison between the truth and the model outputs. Two 10-day forecasts are conducted in order to show the impact of the assimilation on forecasting. The first one takes the model state of the first guess at day 9 as the initial condition and the second one takes that of iteration 13 at day 9 as the initial condition. The SSH RMSE (Fig. 8c) of these two experiments from the truth field shows the improvement of the forecast achieved by improving the forecast initial state using the assimilation procedure. The spatially averaged SSH RMSE reduces from 3.6 cm (the initial guess) to 2.49 cm (the last iteration). 3.3. Experiment three: application to real data The data assimilation technique is applied to a real-world case in this experiment. The model configuration and the initial guess remain the same as in the previous experiments. The observation fields are derived from T/P and Jason-1 along-track SSH data, with data accuracy around 4.2 and 3.9 cm from the two satellite altimeters, respectively. The mapping method is used to interpolate the along-track data to the model grid (Fig. 9a). The mapped altimeter data from 9 January 2004 through 18 January 2004 are assimilated into the NCOM model. The first guess of the model state is chosen arbitrarily and is quite different from the observed field. The variance of the residual SSH field decreases

221

P. Yu et al. / Ocean Modelling 27 (2009) 215–229

SSH "Observations" at Day 0

SSH "Observations" at Day 9

0

24˚N

30˚N

24˚N

2000100

Latitude

30˚N

18˚N

18˚N 98˚W

b

2000100 0

a

90˚W

82˚W

98˚W

90˚W

82˚W

SSH (Iteration 10) at Day 9

SSH (First Guess) at Day 9 30˚N

0

20001000

24˚N 24˚N

2000100

Latitude

30˚N

18˚N 18˚N 98˚W —0.5

—0.3

82˚W

98˚W

—0.1

0.1

SSH RMSE (First Guess)

0.5m

24˚N

18˚N

00

30˚N

20001000

Latitude

0.3

82˚W

SSH RMSE (Iteration 10)

30˚N

24˚N

90˚W

2000 10

c

90˚W

18˚N 98˚W

0.0

90˚W Longitude 0.2

82˚W

98˚W

0.4

0.6

90˚W Longitude 0.8

82˚W

1.0m

Fig. 5. (a) SSH fields from the model control run in idealized experiment one (used as ‘‘observations”) at model day 0 (left) and model day 9 (right); (b) SSH synoptic maps (at model day 9) of the initial guess (left) and iteration 10 (right); (c) The root mean square error (RMSE) of SSH of the initial guess (left) and iteration 10 (right) computed from the control run SSH fields.

to 28% of its initial value following the fifth iteration (Fig. 4c). Visual comparison between the observations and the model outputs (Fig. 9b) and the RMSE map of the model SSH relative to the gridded altimeter data (Fig. 9c) show the improvement of the SSH field by assimilating the observations. The strong anticyclonic eddy in the northwestern GoM, which is in the first guess but not observed by the satellite altimeters, completely disappears and the variability around Cuba is successfully captured. The SSH field inconsis-

tency between model and observations is also evaluated over the along-track observation locations. The along-track SSH RMSE computed along the altimeter ground tracks (Fig. 10) behaves in a manner consistent with the variance of the residual SSH field and drops to 9 cm, 43% of its value from the first iteration (20.7 cm). The improvement is acceptable, given that the adjoint is only based on a reduced-gravity model and the only type of data assimilated are altimetry observations. However, the RMSE is roughly double

222

P. Yu et al. / Ocean Modelling 27 (2009) 215–229

18C Isotherm (”Truth”, Day9)

24˚N

2000 1000

Latitude

30˚N

18˚N 98˚W

90˚W Longitude

82˚W 18C Isotherm (Iter10, Day9)

18C Isotherm (1st Guess, Day9)

2000 1000

Latitude

24˚N

24˚N

2000 1000

30˚N

30˚N

18˚N

18˚N 98˚W

-500

90˚W Longitude -400

82˚W

-300

98˚W

-200

90˚W Longitude -100

82˚W

0m

Fig. 6. The 18 °C isotherm depth (at model day 9) of the ‘‘truth”, or control run (top), the first guess (lower left) and iteration 10 (lower right).

the accuracy of the observations (roughly 4 cm). This limitation in performance is likely due to the representativeness of the subspace constructed from the first baroclinic mode. 4. Discussion It is computationally expensive to implement all but the simplest oceanic data assimilation systems due to the complexity and grid sizes of ocean models, especially high-resolution eddypermitting models. Previous studies have shown success in reducing the dimension of control variables by projecting the state space onto a linear subspace based on statistical decomposition, such as EOF analysis. Due to the much lower degree of freedom in the search space, the computational cost is greatly reduced. However, the efficiency of the assimilation crucially relies on the quality of the representation of the full ocean dynamics by the subspace. The technique developed here starts from a physical point of view and uses dynamical modes as a basis to construct the subspace. The adjoint model is constructed solely based on the first baroclinic mode due to the fact that this mode is dominant in much of the ocean. The vertical normal mode decomposition and reconstruction is used to link the forward and backward models dynamically. The reduced-dynamics variational data assimilation method is applied to the NCOM configured for the domain of the GoM with satellite altimeter data. The system is able to improve the model/ data consistency with a lower computational cost within the 10-

day time assimilation window. For the application to the NCOM simulation used here, the dimension is reduced from 60 vertical layers and five prognostic variables in the forward model to one vertical layer and three variables in the backward reduced-gravity adjoint. The experiments shown here demonstrate that the estimate of the model state, including model SSH and other variables, is improved by applying the new assimilation scheme. A better estimate of the SSH field is transferred throughout the water column through the vertical normal mode decomposition and reconstruction, resulting in an improved estimate of other variables such as the 18 °C isotherm shown in experiment one. However, the system performance is limited by the representativeness of the subspace, e.g., if the first baroclinic mode is not dominant in the region studied, application of this technique may be less efficient or even yields a solution far away from the optimal state. Because of the inconsistency between the forward and backward models, the gradient of the cost function estimated by the adjoint is only approximately correct, which may be the reason for a non-monotonically decreasing cost function in the experiments. This suggests that the convergence may not be guaranteed in some cases, even though the technique has been tested for the difficult case of a first guess of the initial condition very different from that of the control run (idealized experiment one). This inconsistency also induces a limit for the improvement in consistency between the model and observations. A mapping technique has been used for producing regularly gridded SSH data sets from rather sparsely sampled satellite SSH

223

P. Yu et al. / Ocean Modelling 27 (2009) 215–229

SSH "Truth" at Day 0

SSH "Truth" at Day 5 30˚N

30˚N

18˚N

18˚N 82˚W

98˚W

Gridded SSH at Day 0

90˚W

82˚W

24˚N

30˚N

18˚N

98˚W

Gridded SSH — SSH at Day 0

90˚W

82˚W

98˚W

Gridded SSH — SSH at Day 5

18˚N

0

00 100 20

24˚N

18˚N

18˚N 90˚W Longitude

82˚W

—0.5

98˚W

—0.3

82˚W

Gridded SSH — SSH at Day 9

0

0

00 100 20

24˚N

90˚W

30˚N

30˚N

98˚W

20

30˚N

18˚N

98˚W

82˚W

20 10 00 00

24˚N

90˚W

Gridded SSH at Day 9

20 10 00 00

20 10 00 00

Latitude

98˚W

30˚N

18˚N

Latitude

82˚W

Gridded SSH at Day 5

30˚N

24˚N

90˚W

00 100

90˚W

18˚N

20

98˚W

24˚N

24˚N

20

24˚N

00 1 0 0 0

00 1 0 0

0

0

00 1 0 0

24˚N

20

Latitude

30˚N

SSH "Truth" at Day 9

90˚W Longitude —0.1

98˚W

82˚W

0.1

0.3

90˚W Longitude

82˚W

0.5m

Fig. 7. Top row: Control run (‘‘truth”) SSH field from experiment two at Days 0, 5 and 9 of the assimilation time period constructed by applying the CEOF mapping technique to values sampled along locations corresponding to T/P and Jason-1 observation locations (black dots) every ten days. Middle row: SSH fields constructed by applying the mapping technique of Yu et al. (2007) to the model data sampled along synthetic altimeter tracks. Lower row: differences between the model SSH and the mapped synthetic altimeter sampled SSH fields.

observations. Even though this intermediate step is not required, the incorporation of it into the system is beneficial. This technique considers the propagation of information based on the observations from a much longer period compared to the assimilation window and obtains a data set with a finer regular grid. This reduces the duration of the assimilation window, and hence the computational cost, required for the model to capture the propagation of observed features. A shorter assimilation period helps to avoid the occurrence of secondary minima, which may trap the solution away from the primary optimal point (Hoteit, 2008). The experiments are designed to run in a 10-day assimilation window, which is based on the TOPEX/Poseidon ground track repeat cycle. Other timescales over which the method may work have not been determined. Since the mapping technique only considers the first several EOF modes, all the other less important modes are regarded

as noises and ignored, which may help to reduce the noise level of the data assimilated into the models. The forward and backward models are dynamically linked by vertical normal mode decomposition and reconstruction. As a consequence, the assimilation method is able to be designed in a modular fashion. Therefore, it could be applied to other threedimensional ocean models without major changes as long as the first baroclinic mode is dominant in the region being modeled. This feature of the assimilation system is advantageous considering the extreme effort required for developing the adjoint of a complex ocean model. Importantly, when only poor knowledge of the ocean state is known, the approach can be quickly applied to get a better initial condition closer to the truth, at a relatively lower computational cost. The improved initial condition will lead to a better forecast result or can be used as an improved initial guess for more

224

P. Yu et al. / Ocean Modelling 27 (2009) 215–229

a

SSH,Velocity (1st Guess,Day9)

SSH,Velocity (Iter13,Day9)

24˚N

2000

24˚N

10 00

10 00

30˚N

2000

Latitude

30˚N

18˚N

18˚N ••••

98˚W —0.5

b

••••

90˚W —0.3

82˚W

98˚W

—0.1

0.1

0.3

0.5m

30˚N

00 10

24˚N

20

24˚N

20

00 10

00

00

30˚N

18˚N

18˚N 98˚W

90˚W

98˚W

82˚W

SSH RMSE (1st Guess,Forecast)

c

82˚W

SSH RMSE (Iter13,Forecast)

30˚N

30˚N

1 2000 000

24˚N

90˚W

24˚N

18˚N

1 2000 000

Latitude

82˚W

SSH RMSE (Iter13)

SSH RMSE (1st Guess)

Latitude

90˚W

18˚N

98˚W

0.0

90˚W Longitude

82˚W

98˚W

90˚W Longitude

0.2

0.3

0.4

0.1

82˚W

0.5m

Fig. 8. (a) SSH and surface velocity synoptic maps of the initial guess (left) and Iteration 13 (right) at model day 9 from the idealized assimilation experiment two; (b) the RMSE of SSH of the initial guess (left) and iteration 13 (right) computed from the control run SSH fields; (c) the RMSE of SSH following a 10-day forecast initialized from the first guess (left) and from iteration 13 (right).

complicated assimilation procedures to increase the speed of their convergence. The main goal of this study is to provide a new method for estimating the ocean state with a lower computational cost compared to that usually required when solving complicated 4D-Var problems. The method has been applied to a model of the GoM that has high resolution in all three spatial dimensions. Experiments have been conducted to test the effectiveness of the method, and to demonstrate its applicability to a real-world case by assimilating an SSH data set produced from satellite altimeter data using a new

gridding technique. Limitations of this data assimilation technique, and assumptions that have necessarily been made to apply the method include: the representation of full model dynamics by a first baroclinic mode approximation for the adjoint, the restriction that the method only be applied where the SSH variability is dominated by the first baroclinic mode, the assumption of a flat bottom for the computation of the vertical modal structure, and the restriction that only SSH data be assimilated. Additionally, particular choices for processing of the SSH data and mapping density onto temperature and salinity profiles have been made

225

P. Yu et al. / Ocean Modelling 27 (2009) 215–229

a

Gridded SSH (Day0)

Gridded SSH (Day9)

24°N

18°N

20

24°N

00 10 00

30°N

200 100 0 0

Latitude

30°N

18°N

98°W

b

90°W

82°W

98°W

30°N

00 10

00 10 2000

24°N

2000

Latitude

30°N

18°N

18°N

1m/s

1m/s

98°W —1.0

c

90°W

82°W

—0.6

98°W

—0.2

0.2

SSH RMSE (1st Guess)

0.6

82°W 1.0m

24°N

18°N

200 10 0 00

30°N

200 10 0 00

Latitude

90°W

SSH RMSE (Iter14)

30°N

24°N

82°W

SSH,Velocity (Iter14,Day9)

SSH,Velocity (1st Guess,Day9)

24°N

90°W

18°N

98°W

0.0

90°W Longitude

82°W

0.2

0.4

98°W

0.6

90°W Longitude 0.8

82°W

1.0m

Fig. 9. (a) Synoptic maps of the gridded SSH fields produced applying the CEOF mapping method of Yu et al. (2007) to T/P and Jason-1 data shown at model day 0 (January 9, 2004) and model day 9 (January 18, 2004); (b) synoptic maps of SSH and surface velocity for the initial guess (left) and iteration 14 (right) at model day 9; (c) the RMSE of SSH of the initial guess (left) and iteration 14 (right) computed with respect to the gridded observation fields.

for implementation of the method in this work. However, the performance of the assimilation system can likely be improved by modification of various algorithms used in the application of this technique with other options, further highlighting the advantage of the modularity of the method.

chael Navon and Xiaolei Zou for their advice. Funding for this project was provided by the Office of Naval Research, the NOAA Office of Global Programs, and the NASA Physical Oceanography Program.

Acknowledgements

Appendix A. Vertical normal mode decomposition from the NCOM governing equations

Suggestions and comments from the two anonymous reviewers are greatly appreciated. The authors thank Drs. Paul Martin and Alan Wallcraft for assistance with the NCOM development, and Drs. I. Mi-

Following Philander (1990), a method of computing vertical normal modes is derived from the NCOM model equations. The NCOM equations in continuous form are given by

226

P. Yu et al. / Ocean Modelling 27 (2009) 215–229

By again using equation (A8) it can be shown that

0.3

1 0

q rh ðAH rh qÞ ¼ arh ðAH rh TÞ þ brh ðAH rh SÞ

0.25

ðA12Þ

and

RMSE

0.2

q1 0 ðK H qz Þz ¼ aðK H T z Þz þ bðK H Sz Þz

0.15

ðA13Þ

Substituting (A12) and (A13) into (A11) gives a prognostic equation for density

0.1

1 qt =q0 ¼ q1 0 r  ðvqÞ þ QðaT þ bSÞ þ q0 rh ðAH rh qÞ þ q1 0 ðK H qz Þz  aQ r cz

0.05 0 1

2

3

4

5

6

7 8 9 10 11 12 13 14 15 Iteration

Fig. 10. The root mean square error (RMSE) of model SSH anomaly from the alongtrack satellite data computed along the satellite ground tracks shown.

  @u 1 @p @ @u ¼ r  ðvuÞ þ Qu þ fv  þ Fu þ KM ; @t @z @z q0 @x   @v 1 @p @ @v ; ¼ r  ðvvÞ þ Qv  fu  þF þ K @t q0 @y v @z M @z @p ¼ qg; @z @u @v @w þ þ ¼ Q; @x @y @z   @T @ @T @c ¼ r  ðvTÞ þ QT þ rh ðAH rh TÞ þ KH þ Qr ; @t @z @z @z   @S @ @S ¼ r  ðvSÞ þ QS þ rh ðAH rh SÞ þ KH ; @t @z @z

ðA1Þ ðA2Þ ðA3Þ ðA4Þ ðA5Þ ðA6Þ

q ¼ qðT; S; zÞ;

ðA7Þ

where t is the time, x, y, and z are the three coordinate directions, u, v, and w are the three components of the velocity, Q is a volume flux source term, v is the vector velocity, T is the potential temperature, S is the salinity, rh is the horizontal gradient operator, f is the Coriolis parameter, p is the pressure, q is the water density, q0 is a reference water density, g is the acceleration of gravity, F u and F v are the horizontal mixing terms for momentum, AH is the horizontal mixing coefficient for scalar fields (potential temperature and salinity), K M and K H are the vertical eddy coefficients for momentum and scalar fields, respectively, Q r is the solar radiation, and c is a function describing the solar extinction. The model equations are simplified as follows. First, the equation of state is simplified to a linear form as

q ¼ q0 ½1  aðT  T 0 Þ þ bðS  S0 Þ;  q1 @@Tq, 0

ðA8Þ

1 @q q0 @S

where a ¼ and b ¼ . So the time dependence of density can be written as

ðA9Þ

  @u 1 @p @ A @u ¼ fv  þ @t q0 @x @z N2 @z   @v 1 @p @ A @v ¼ fu  þ @t q0 @y @z N2 @z @p ¼ qg @z @u @v @w þ þ ¼0 @x @y @z @q ¼ r  ðvqÞ þ rh ðAH rh qÞ @t   @ A @q @q  w þ @z N2 @z @z   @ A @q þ ; i:e:; @z N2 @z   @q N2 @ A @ q ¼ q0 w þ 2 @z N @z @t g

      1 @q @T @S þ r  ðvTÞ þ b þ r  ðvSÞ þ r  ðvqÞ ¼ a @t @t q0 @t ðA10Þ allows equations (A5) and (A6) to be substituted to obtain the form

  @q þ r  ðvqÞ ¼ a½QT þ rh ðAH rh TÞ þ ðK H T z Þz þ Q r cz  q0 @t 1

where the subscript z stands for @[email protected] and the subscript t stands for @[email protected]

ðA16Þ ðA17Þ ðA18Þ

ðA19Þ

  @ A p @z N 2 zz

ðA20Þ

The method of separation of variables is applied by assuming M X

U m ðx; y; tÞRm ðzÞ

ðA21Þ

V m ðx; y; tÞRm ðzÞ

ðA22Þ

m¼0



M X m¼0

XM p ¼ g ðx; y; tÞRm ðzÞ m¼0 m g q0 w¼

ðA11Þ

ðA15Þ

q is eliminated by substituting (A17) into (A19) yielding



Expanding the material time derivative as

þ b½QS þ rh ðAH rh SÞ þ ðK H Sz Þz ;

that can replace Eqs. (A5)–(A7). Eq. (A14) and Eqs. (A1)–(A4) are the reduced set of equations. The problem needs to be further simplified by linearization of the above equations. To linearize these equations it is necessary to assume that the ocean, in its basic state, is motionless and has  , which is a function of depth only. Motion and its density q ¼ q associated density variation q0 ðx; y; z; tÞ are considered as small perturbations to this basic state. For the sake of mathematical convenience, it is assumed that K M ¼ K H ¼ A=N 2 , where N is de z =q0 Þ1=2 (the Brunt-Väisälä frequency) and A is fined as N ¼ ðg q a constant. With these assumptions and neglecting the source term Q and solar radiation Qr, the linearized unforced NCOM equations become

pzt þ q0 wN2 ¼

  dq @ q dT @ q dS dT dS þ  ¼ aq0 þ bq0 : ¼  dt dt dt @T S dt @T T dt

ðA14Þ

M X

W m ðx; y; tÞSm ðzÞ

ðA23Þ ðA24Þ

m¼0

where m is the mode number and M is the upper limit for m which goes to infinity for continuous functions.

227

P. Yu et al. / Ocean Modelling 27 (2009) 215–229

with boundary conditions

Substitution yields

ðU t  fV þ g gx ÞR ¼ ðV t þ fU þ g gy ÞR ¼



A



URz z VRz 2

N2  A

N ðU x þ V y ÞR þ WSz ¼ 0   WN2 A gt Rz þ g R S¼ zz g N2 z

ðA25Þ ðA26Þ

z

ðA27Þ ðA28Þ

where the subscript m has been dropped. By rewriting (A27) as

U x þ V y Sz 1 ¼ ¼ ; gh W R

R ¼ ghSz

ðA30Þ

Separating the z-dependence in (A25) and (A26) yields

N2

¼ z

R gh

Sz þ gh

Szz

N

S þ gh

N

2

z

ðA33Þ

¼0

Z

0

Rm Rn dz ¼ 0

m–n

ðA48Þ

The baroclinic modes can be obtained as 0

uRn dz ¼

Z

0

"

X

H

# U m Rm Rn dz ¼

m

X

Um

Z

m

0

Rm Rn dz ¼ U n

H

Z

0

H

R2n dz

ðA49Þ Z

Un ¼

Z uRn dz

0 H

0

H

R2n dz

ðA50Þ

R2n dz

ðA51Þ

Similarly,

ðA34Þ

z

¼C

ðA47Þ

Based on the Sturm–Liouville general solution, the eigenfunctions, Rn, form a complete set and are orthogonal in the sense that



2

R ¼0 gh

So,

which is integrate vertically to arrive at the relation

Szz

þ

N2

ðA32Þ

Plugging (A30) into (A33) yields





Rz

H

Both (A31) and (A32) lead to

Rz



Z

2



ðA46Þ

Eq. (A43) along with the boundary conditions (A45) and (A46) is a Sturm–Liouville problem, and the discrete version is an eigenvalue problem, which can be solved via standard numerical methods. (A40)–(A42) are governing equations for each mode, which are shallow-water equations. From (A43) and (A44), the equation for R is

ðA31Þ

and

V t þ fU þ g gy ðRz =N Þz 1 ¼ ¼ gh AV R

ðA45Þ

z¼0

H

2

U t  fV þ g gx ðRz =N Þz 1 ¼ ¼ gh AU R

z ¼ zbottom

S  hSz ¼ 0;

ðA29Þ

where h is the constant of separation, one gets



St ¼ 0;

Vn ¼

gn ¼

Z

0

vRn dz

H Z 0

ðA35Þ

0

H

p g q0

H

Z

Rn dz

Z

0

H

R2n dz

ðA52Þ

where C is an arbitrary constant. Assuming S ¼ 0 and Szz ¼ 0 at the bottom leads to C ¼ 0 and (A35) becomes

Appendix B. Derivation of the adjoint of the reduced-gravity model

Szz

Adjoint equations are derived from the equations for a reducedgravity ocean model, using the calculus of variations and forming a Lagrange cost function by adding the constraints multiplied by the Lagrange multipliers. The cost function is defined as

N

2

þ

S ¼0 gh

ðA36Þ

The boundary conditions are

w¼S¼0 z ¼ zbottom pt w ¼ 0; z ¼ 0: q0 g

ðA37Þ ðA38Þ

F ¼ Km

Substituting (A23) and (A24) into (A38) leads to

Z

 ðA40Þ þ

kv

a

  ðA41Þ þ kg  ðA42Þ dxdydt:

ðB1Þ

Assuming the observation errors at different locations are uncorrelated and constant observation error variances reduces Km to a constant. So (B1) can be rewritten as

So,

ðA39Þ

Thus, the vertical normal mode decomposition problem is reduced to the following equations:

U t  fV þ g gx ¼ AU=gh V t þ fU þ g gy ¼ AV=gh

ðA40Þ ðA41Þ

gt þ hðU x þ V y Þ ¼ Ag=gh

ðA42Þ

N2 Szz þ S¼0 gh R ¼ ghSz

ku

a

x;y;t

WSz ¼ 0: Ux þ V y

S  hSz ¼ 0:

ðgðx; y; tÞ  gobs ðx; y; tÞÞ2 dxdydt



Z

þ

WS  gt ¼ 0 i:e:;WS  gt

x;y;t

ðA43Þ ðA44Þ



1

Z

a þ þ

x;y;t

1

ku ½U t  fV þ g gx þ AU=ghdxdydt

Z

a Z

þ Km

kv ½V t þ fU þ g gy þ AV=ghdxdydt

ðIIÞ

kh ½gt þ hðU x þ V y Þ þ Ag=ghdxdydt

ðIIIÞ

ðgðx; y; tÞ  gobs ðx; y; tÞÞ2 dxdydt

ðIVÞ

x;y;t

x;y;t

ðIÞ

Z x;y;t

ðB2Þ

where ku , kv , and kh are the adjoint variables with respect to the control variables, U, V, and g. At the minimum of the cost function,

228

P. Yu et al. / Ocean Modelling 27 (2009) 215–229

F has a stationary point, and its first derivative with respect to all control variables must vanish. In order to find the adjoint equations, one has to set the first derivative of the associated Lagrange function zero, i.e.,

So, the governing equation for kv is

@F ¼ 0; @U

ðB3Þ

Consider (B5),

@F ¼ 0; @V

ðB4Þ

1 @kv

a @t @ 1 @g a

@F ¼ 0; @g

ku ½g gx dxdydt

x;y;t

Z  Z 1 @ @ @ku g ðku gÞdxdydt  gg dxdydt a @ g x;y;t @x @x x;y;t Z 1 @ku ¼ g dxdydt a x;y;t @x

ðB5Þ

(B3) is considered first:

Z 1 @ @ ku ½U t dxdydt ¼ ðk UÞdxdydtt a @U x;y;t @t u x;y;t  Z @ku 1 @ku U dxdydt ¼  dxdydt @t a x;y;t @t x;y;t

@ 1 @U a Z 

@ 1 @U a @ 1 @U a



Z ku x;y;t

Z

 Z AU 1 Aku dxdydt ¼ dxdydt gh a x;y;t gh

kv ½fUdxdydt ¼

x;y;t

1

Z

a

@ 1 @g a

@ @g

ðB8Þ

Z

@ @g ðB9Þ

1 @ku 1 @kh 1 Aku  fk þ h ¼ a @t a v @x a gh

ðB10Þ

ku ½fVdxdydt ¼ 

x;y;t

@ 1 @V a

Z

f ku dxdydt

ðB11Þ

x;y;t

Z

¼

1

x;y;t

Z

a

x;y;t



Z kv x;y;t

Z

@ ðkv VÞdxdydt  @t

Z x;y;t

V

@kv dxdydt @t



@kv dxdydt @t 

AV 1 dxdydt ¼ gh a

x;y;t

Z kh x;y;t

@ Km @g

  Z Ag Akh dxdydt ¼ dxdydt gh x;y;t gh

Z x;y;t

½g  gobs 2 dxdydt ¼ 2K m

Z x;y;t

ðg  gobs Þdxdydt

ðB18Þ

ðB19Þ

ðB20Þ

  @kh g @ku @kv Akh þ þ ¼ þ 2K m ðg  gobs Þ @t a @x @y gh

ðB21Þ

Akv dxdydt gh

@ku @kh Aku  f kv þ ah ¼ @t @x gh

ðB22Þ ðB23Þ ðB24Þ

For the scaling purposes, let a ¼ g=h, and the equations become

ðB12Þ Z

x;y;t

@kv @kh Akv þ f ku þ ah ¼ @t @y gh   @kh g @ku @kv Akh þ þ ¼ þ 2K m ðg  gobs Þ @t a @x @y gh

x;y;t

a @V

@ 1 @V a

a

Z

kv ½V t dxdydt

1 @

¼

1

Z

Reorganize the above equations

Similarly, when (B4) is considered,

Z

ðB17Þ

So, the governing equation for kh is

So, the governing equation for ku is

@ 1 @V a

x;y;t

kh ½gt dxdydt Z  Z @ @k ðkh gÞt dxdydt  g h dxdydt ¼ @ g x;y;t @t x;y;t Z @kh ¼ dxdydt: x;y;t @t

x;y;t

kh ½hU x dxdydt x;y;t  Z  Z @ @kh ðkh UÞx dxdydt  h U ¼ h dxdydt @U @x x;y;t x;y;t Z @kh ¼ h dxdydt: @x x;y;t

Z

kv ½g gy dxdydt Z  Z 1 @ @ @kv g ðkv gÞdxdydt  gg ¼ dxdydt a @ g x;y;t @y @y x;y;t Z 1 @kv ¼ g dxdydt a x;y;t @y

ðB6Þ

ðB7Þ

f kv dxdydt

ðB16Þ

Similarly,

Z

where the initial conditions, ku jt¼s ¼ 0 and Ujt¼0 are used

@ @V

Z

ðB15Þ

¼

and

@ @U

1 @kv 1 Akv þ f ku þ h ¼ a @y a gh

ðB13Þ

@ku @kh Aku  f kv þ g ¼ @t @x gh @kv @kh Akv þ f ku þ g ¼ @t @y gh   @kh @ku @kv Akh þh þ ¼ þ 2K m ðg  gobs Þ @t @x @y gh

ðB25Þ ðB26Þ ðB27Þ

kh ½hV y dxdydt

x;y;t

 Z  Z @ @kh ðkh VÞy dxdydt  h V h dxdydt ¼ @V @y x;y;t x;y;t Z @kh h dxdydt ¼ @y x;y;t

References

ðB14Þ

Barnett, T.P., 1983. Interaction of the monsoon and Pacific trade wind system at interannual time scales. Part I: the equatorial zone. Mon. Weather Rev. 111, 756–773. Bennett, A.F., 2002. Inverse Modeling of the Ocean and Atmosphere. Cambridge University Press. 234 pp.

P. Yu et al. / Ocean Modelling 27 (2009) 215–229 Blayo, E., Dubiano, S., Vidard, P.A., Le Dimet, F.X., 2003. Reduced order strategies for variational data assimilation in oceanic models. In: Sportisse, B., Le Dimet, F.-X. (Eds.), Data Assimilation for Geophysical Flows. Springer-Verlag. Chua, B., Bennett, A.F., 2001. An inverse ocean modeling system. Ocean Modell. 3, 137–165. Cooper, M., Haines, K., 1996. Altimetric assimilation with water property conservation. J. Geophys. Res. 101, 1059–1077. Courtier, P., Talagrand, O., 1990. Variational assimilation of meteorological observations with the direct and adjoint shallow-water equation. Tellus 42A, 531–549. DaSilva, A., Young, A.C., Levitus, S., 1994. Atlas of Surface Marine Data 1994, Algorithms and Procedures, NOAA Atlas NESDIS 6, Natl. Oceanic Atmos. Admin., vol. 1, Silver Spring, Md. Dewar, W.K., Morris, M.Y., 2000. On the propagation of baroclinic waves in the general circulation. J. Phys. Oceanogr. 30, 2637–2649. Dewar, W.K., Huang, R.X., 2001. Adjustment of the ventilated thermocline. J. Phys. Oceanogr. 31, 1676–1697. Evensen, G., 1994. Sequential data assimilation with a nonlinear quasi-geostropic model using Monte Carlo methods to forecast error statistics. J. Geophys. Res. 105, 19477–19498. Fukumori, I., 2002. A partitioned Kalman filter and smoother. Mon. Weather Rev. 130, 1370–1383. Fukumori, I., Malanotte-Rizzoli,, P., 1995. An approximate Kalman filter for ocean data assimilation: an example with an idealized Gulf Stream model. J. Geophys. Res. 100, 6777–6793. Gill, M., Malanotte-Rizzoli, P., 1991. Data assimilation in meteorology and oceanograph. Adv. Geophys. 33, 141–266. Ghil, M., Cohn, S.E., Tavantzis, J., Bube, K., Isaacson, E., 1981. Application of estimation theory to numerical weather prediction. In: Bengtsson, L., Ghil, M., Kallen, E. (Eds.), Dynamical Meteorological: Data Assimilation Methods. Springer-Verlag, pp. 139–224. Hoteit, I., 2008. A reduced-order simulated annealing approach for fourdimensional variational data assimilation in meteorology and oceanography. Int. J. Numer. Meth. Fluids. doi:10.1002/fld.1794. Hoteit, I., Cornuelle, B., Kohl, A., Stammer, D., 2005. Treating strong adjoint sensitivities in tropical eddy-permitting variational data assimilation. Q.J.R. Meteorol. Soc. 131, 3659–3682. Kohl, A., Willebrand, J., 2002. An adjoint method for the assimilation of statistical characteristics into eddy-resolving ocean models. Tellus 54A, 406–425. Le Dimet, F.X., Navon, I.M. 1988. Variational and optimization methods in meteorology: a review. Technical Report FSU-SCRI-88-144, Florida State University, Tallahassee, Florida, 83 pp. Le Dimet, F.X., Talagrand, O., 1986. Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects. Tellus 38A, 97–110. Liu, D.C., Norcedal, J., 1989. On the limited memory BFGS method for large scale optimization. Math. Prog. 45, 503–528.

229

Liu, Z., 1999. Planetary wave modes in thermocline: non-Doppler-shift mode, advective mode and Green mode. Q.J.R. Meteorol. Soc. 125, 1315–1339. Martin, P.J., 2000. A description of the Navy Coastal Ocean Model Version 1.0. NRL/ FR/7322-009962, Naval Research Laboratory, Stennis Space Center, Mississippi, 39 pp. Morey, S.L., Schroeder, W.W., O’Brien, J.J., Zavala-Hidalgo, J., 2003. The annual cycle of riverine influence in the eastern Gulf of Mexico basin. Geophys. Res. Lett. 30, 1867–1870. Morey, S.L., J. Zavala-Hidalgo, J.J. O’Brien, 2005. The seasonal variability of continental shelf circulation in the northern and western Gulf of Mexico from a high-resolution numerical model in Circulation of the Gulf of Mexico: observations and Models, In: Sturges, W., Lugo-Fernandez, A. (Eds.), Geophys. Mongr. Ser. 161, AGU, Washington, DC. doi: 10.1029/161GM16. Navon, I.M., Daescu, D.N., Liu, Z., 2005. The impact of background error on incomplete observations for 4D-Var data assimilation with the FSU GSM. In: Sunderam V.S. et al. (Eds.), Computational Science-ICCS 2005, LNCS, vol. 3515, Springer-Verlag, pp. 837–844. Navon, I.M., Legler, D.M., 1987. Conjugate-gradient methods for large-scale minimization in meteorology. Mon. Weather Rev. 115, 1479–1502. Navon, I.M., Zou, X., Derber, J., Sela, J., 1992. Variational data assimilation with an adiabatic version of the NMC spectral model. Mon. Weather Rev. 122, 1433– 1446. ˇ a, and the Southern Oscillation. Academic Press, Philander, S.G., 1990. El Ninˇo, La Nin Inc., pp. 160–167. Robert, C., Dubiano, S., Blayo, E., Verron, J., Blum, J., Le Dimet, F.X., 2005. A reducedorder strategy for 4D-Var data assimilation. J. Mar. Syst. 57, 70–82. Siegel, D.A., Mcgillicuddy Jr., D.J., Fields, E.A., 1999. Mesoscale eddies, satellite altimetry, and new production in the Sargasso Sea. J. Geophys. Res. 104, 13359– 13379. Shriver, J.F., Johnson, M.A., O’Brien, J.J., 1991. Analysis of remotely forced oceanic Rossby waves off California. J. Geophys. Res. 96, 749–757. Shu, L., Clarke, A.J., 2002. Using an ocean model to examine ENSO dynamics. J. Phys. Oceanogr. 32, 903–923. Wunsch, C., 1996. The Ocean Circulation Inverse Problem. Cambridge University Press, 437 pp. Wunsch, C., 1997. The vertical Partition of oceanic horizontal kinetic energy. J. Phys. Oceangr. 27, 1770–1794. Wunsch, C., Heimbach, P., 2007. Practical global oceanic state estimation. Physica D 230, 197–208. Yu, L., O’Brien, J.J., 1991. Variational estimation of the wind stress drag coefficient and the oceanic eddy viscosity profile. J. Phys. Oceanogr. 21, 709–719. Yu, P., Morey, S.L., Zavala-Hidalgo, J., 2004. New mapping method to observe propagating features. Sea Technology 45 (5), 20–24. Zou, X., Navon, I.M., Berger, M., Phua, K.H., Schlick, T., Le Dimet, F.X., 1993. Numerical experience with limited-memory quasi-Newton and truncated Newton methods. SIAM J. Numer Optimization 3, 582–608.