Application of 4D-Variational data assimilation to the California Current System

Dynamics of Atmospheres and Oceans 48 (2009) 69–92

Contents lists available at ScienceDirect

Dynamics of Atmospheres and Oceans journal homepage: www.elsevier.com/locate/dynatmoce

Application of 4D-Variational data assimilation to the California Current System G. Broquet a,∗, C.A. Edwards b, A.M. Moore b, B.S. Powell a, M. Veneziani a, J.D. Doyle c a b c

University of California Santa Cruz, Institute of Marine Sciences, 1156 High Street, Santa Cruz, CA 95064, United States University of California Santa Cruz, Department of Ocean Sciences, CA 95064, United States Naval Research Laboratory, 7 Grace Hopper Ave., Stop 2 Monterey, CA 93943-5502, United States

a r t i c l e

i n f o

Article history: Available online 5 April 2009

Keywords: In situ/surface data assimilation 4DVAR California Current System ROMS Hindcasts/analysis

a b s t r a c t The Incremental Strong constraint 4D-Variational (IS4DVAR) data assimilation system of the Regional Ocean Model System (ROMS) is used to study the controllability of a realistic, high resolution configuration of the California Current System. The configuration and results of assimilating both satellite-derived surface observations along with in situ data are presented. Results show consequent improvements in many characteristics of the model circulation, and some of the strengths of the adjoint method for data assimilation are highlighted. General issues of the sensitivity of the results to the configuration of ROMS-IS4DVAR are also discussed. © 2009 Elsevier B.V. All rights reserved.

1. Introduction The oceanic circulation along the U.S. west coast is dominated by the California Current System (CCS; Hickey, 1979), a set of currents along the oceanic eastern boundary, consisting of the broad equatorward California Current (CC), the poleward California Under Current (CUC), the poleward Davidson Current (DC), and the coastal alongshore jet associated with periods of upwelling and downwelling. This system has been intensively observed in the upper 1000 m, particularly through the California Cooperative Ocean Fisheries Investigation (CalCOFI; Chereskin and Trunnell, 1996) leading to a good description of some of the major circulation features. Theoretical and numerical modeling have led to significant

∗ Corresponding author. Tel.: +1 831 459 1306; fax: +1 831 459 4882. E-mail address: [email protected] (G. Broquet). 0377-0265/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.dynatmoce.2009.03.001

70

G. Broquet et al. / Dynamics of Atmospheres and Oceans 48 (2009) 69–92

progress in the understanding of dynamical mechanisms of the CCS (McCreary et al., 1987; Batteen, 1997; Hickey, 1998). By definition, data assimilation is the optimal way of combining both observations and models. Data assimilation methods have reached a level of maturity that warrants their use towards an improved knowledge of the unresolved dynamics in the CCS and towards operational forecasting. The present study relies on the association of a recently developed model and data assimilation system that have already shown great promise, and that constitute a new approach to the study of the CCS. Through the development of nesting techniques and coastal modeling, recent efforts for modeling the CCS in realistic configurations have focused on local areas of the U.S. west coast, such as the northern California coast with the Princeton Ocean Model (POM) (Gan and Allen, 2002a, b), the Monterey Bay area (Shulman et al., 2007) with the Navy Coastal Ocean Model (NCOM), using the Regional Oceanic Modeling System (ROMS) in the Southern California Bight (SCB; Di Lorenzo, 2003), and using ROMS in the northern California (Cervantes and Allen, 2006). In accordance with observations, these models have revealed the primary influence of atmospheric wind forcing, bottom topography, and the shape of the coastline on the dynamics and ocean variability. The ROMS configurations of Marchesiello et al. (2003), covering the entire U.S. west coast and CCS, demonstrated the importance of the variability forced by internal dynamics such as the interactions with the bathymetry and baroclinic instabilities. This work also highlighted the skill of ROMS for regional modeling, especially near the continental shelf (Haidvogel et al., 2008). Veneziani et al. (2009a) also used ROMS to develop a CCS configuration, but using realistic rather than climatological data to force the model at the surface and at the open boundaries. These authors showed that realistic forcing substantially improves both the mean circulation and the eddy kinetic energy variability near the coast, especially near capes. Data assimilation experiments have also been conducted in some local areas of the CCS. Sequential filter/Optimal Interpolation (OI) systems based on the Kalman Filter (KF) have been used in a midOregon shelf configuration to assimilate moored velocity profiles (Kurapov et al., 2005a, b) and surface current estimates from HF radars (Oke et al., 2002). In the Monterey Bay area, a system also derived from the sequential scheme of the KF via the Error Subspace Statistical Estimation (ESSE) by Lermusiaux (2007), has been tested for real-time predictions. Nudging techniques have also been applied in the Monterey bay area (Shulman et al., 2002). The recent application of a 3D-Variational (3DVAR) system to the configuration of Marchesiello et al. (2003)(Li et al., 2008a,b) is the first attempt to use data assimilation for controlling the whole CCS. The ROMS-4DVAR system has been developed recently using extensions of ROMS, such as the tangent linear model (TLROMS) of the original non-linear ROMS model (NLROMS), and the adjoint (ADROMS) of TLROMS, which are fundamental tools that provide opportunities for conducting various types of analysis on complex ocean circulations (Moore et al., 2004). The ROMS-4DVAR data assimilation system has already been evaluated in an idealized coastal upwelling environment that captures some of the important dynamics of the CCS (Di Lorenzo et al., 2007). The Incremental Strong constraint version of ROMS-4DVAR (ROMS-IS4DVAR) has been applied in an operational real-time forecasting system in the Intra-Americas Sea by Powell et al. (2008). Sensitivity studies, using ADROMS, a natural prerequisite for any application of ROMS-IS4DVAR, were conducted for the SCB by Moore et al. (2009) and for the CCS by Veneziani et al. (2009b), demonstrating the robustness of the ROMS adjoint-based tools, and the possibility of controlling the model using the ROMS-IS4DVAR system. The CCS configuration of Veneziani et al. (2009a) was mainly developed to explore the sensitivity of the regional circulation to various surface atmospheric forcing and open boundary conditions. Application of ROMS-IS4DVAR with its capability to adjust the initial conditions, surface forcing and open boundary conditions, will provide a very well adapted tool for these studies. The application of the system with the adjustment of initial conditions only is a first step toward exploring the possibility of controlling the circulation using a limited set of observations. The objective of this paper is to present the main aspects of the application of ROMS-IS4DVAR to the CCS and the overall improvements in the circulation that result from this application. In the following, we will describe in Section 2 the models of the CCS used, and present in Section 3 an overview of the theory underpinning the ROMS-IS4DVAR system. The set up of the ROMS-IS4DVAR system as applied to the CCS and the related experimental protocol is described in Section 4. The results obtained from

G. Broquet et al. / Dynamics of Atmospheres and Oceans 48 (2009) 69–92

71

Fig. 1. Bathymetry of (a) WC3 and (b) WC10 in m depth.

the experiments and the system behavior are described in Section 5. Conclusions are given in Section 6. 2. Model and configuration ROMS is a free-surface, hydrostatic, primitive equation model using a mode-splitting treatment (between the barotropic and baroclinic modes) for the time-stepping. It uses terrain-following vertical coordinates, which is advantageous for regional applications, and is particularly adapted for high resolution coastal zones (Shchepetkin and McWilliams, 2005; Haidvogel et al., 2008). However, the use of this vertical coordinate is predicated on smoothing of the bathymetry which can affect the coastal circulation. Particular attention is paid to minimizing the known deficiencies and errors in the horizontal pressure-gradient (Shchepetkin and McWilliams, 2003). A large and diversified community of users has applied ROMS successfully to model both coastal and basin-scale circulations.1 Sub-grid scale parameterizations of vertical mixing in ROMS can be modeled using different options such as the k–kl (Mellor–Yamada Level 2.5), the k–  and the k– ω turbulence closures represented as parameterizations of the Generic Length Scale (GLS) vertical mixing scheme (Warner et al., 2005), and the KPP mixing scheme of Large et al. (1994). Several of these options are used here because vertical diffusion in NLROMS determines the diffusion parameters used by TLROMS and ADROMS, and this choice appears to significantly influence the vertical data extrapolation during data assimilation. The ROMS configurations used for the CCS cover the domain 134–115.5◦ W and 30–48◦ N with Mercator grids of resolution 1/3◦ and 1/10◦ , and respectively 30 and 42 vertical levels (this yields a vertical resolution of ∼0.3–8 m along the continental shelf and ∼300 m in the deep ocean offshore). The 1/3◦ model was used for studying the effects of different configurations of ROMS-IS4DVAR on the resulting circulations. The two configurations will be respectively referred to as West Coast 1/3◦ (WC3) and West Coast 1/10◦ (WC10). WC10 is described in detail by Veneziani et al. (2009a,b). Fig. 1 shows the domains and bathymetry for WC3 and WC10, which were both generated using ETOPO2 and smoothed with a Shapiro filter, to remove overly large gradients that are problematic for the model. The bathymetry is characterized by the continental shelf, narrow in the original data but broadened with the application of the Shapiro filter, bordering a coastline dominated by several capes that are not sufficiently resolved by WC3. The model was run in both configurations using harmonic horizontal viscosity (mixing coefficient: 4 m2 s−1 ), free-slip boundary conditions, a quadratic bottom friction (bottom drag coefficient: 2.5 × 10−3 ) and alternately k- GLS, k-ω GLS, and KPP vertical mixing schemes.

1 Some of the major applications using the three main versions of ROMS are described at http://www. myroms.org/papers, http://www.atmos.ucla.edu/cesr/ROMS page.html#Selected References and http://www.brest.ird.fr/ Roms tools/references.html.

72

G. Broquet et al. / Dynamics of Atmospheres and Oceans 48 (2009) 69–92

Fig. 2. Summer (July–August–September) mean for 2000–2004 of meridional velocity in m s−1 —vertical cross-shore section at latitude 37◦ N in (a) FREE3 and (b) FREE10.

Open boundary conditions were imposed using monthly averaged data from the Estimating the Circulation and Climate of the Ocean-Global Ocean Data Assimilation Experiment (ECCO-GODAE) model (Heimbach et al., 2006), with a 1◦ resolution horizontally and 23 vertical geopotential z levels vertically. The formulation of the open boundary conditions (OBC) in ROMS is described in Marchesiello et al. (2001). For the barotropic mode, the Flather (1976) condition completed with an Orlanski radiation condition (to estimate separately the free-surface height) as proposed in Chapman (1985) is used. For the baroclinic mode, practical aspects of TLROMS and ADROMS calculations prevent us from using the radiation OBC. Instead, clamped boundary conditions are used along with a sponge layer (with higher viscosity) near the boundary to attenuate any wave reflections produced by the clamped conditions. River forcing was not included in the boundary conditions. The surface forcing data are provided daily by the multi-grid atmospheric component of the Coupled Ocean-Atmosphere Mesoscale Prediction System (COAMPS) model (Hodur et al., 2002; Doyle et al., 2009). COAMPS has a 3–9 km resolution along the California and Oregon coasts. Air temperature, surface pressure, humidity, wind speed and direction, short- and long-wave radiation, and precipitation from COAMPS are used to compute the ROMS surface boundary conditions with the bulk flux formulation of Fairall et al. (1996). The data assimilation experiments are based on NLROMS simulations that are initilized on January 1, 1999 after 7 years of spin-up with climatological forcings, Levitus climatology for OBC and Comprehensive Ocean-Atmosphere Data Set (COADS) data for the surface forcings (the spin-up is initialized with the interpolation of a snapshot from the ECCO-GODAE model). The NLROMS simulation on WC3 (WC10) will be referred to as FREE3 (FREE10). Both FREE3 and FREE10 were run for the 6-year period 1999–2004. Fig. 2 shows vertical sections of July–September mean meridional velocities along 37◦ N close to Monterey Bay, and indicates that both FREE3 and FREE10 are able to reproduce many mean features of the CCS. In both cases the equatorward CC in the upper few hundreds meters, the strong northward CUC near the coast topography, and the nearshore equatorward coastal jet reaching 100 m depth are

G. Broquet et al. / Dynamics of Atmospheres and Oceans 48 (2009) 69–92

73

Fig. 3. (a and b) June SST mean in ◦ C and (c and d) summer SSH mean in m for the period 2000–2004 in (a and c) FREE3 and (b and d) FREE10.

captured. A strong poleward deep and broad subsurface current abutting the CUC, identified in other studies such as in Marchesiello et al. (2003), is present in both configurations. Important differences are noticeable between the two simulations, such as the position of the CUC core and the mean intensity and depth of the coastal jet. Fig. 3 shows the mean SST for June and the July–September mean SSH for FREE3 and FREE10. FREE3 reproduces the main mean features of the CCS and the seasonal variability, with upwelling occurring for example with comparable signatures in both configurations. Meanders of the coastal jet can even be seen in WC3 even though the major coastal capes are not well resolved. Discontinuities between the internal circulation and the clamped boundary conditions are more visible in FREE3 than in FREE10 because they occur at the larger scale of the WC3 resolution. The dome of spurious SSH in the vicinity of the northern boundary (Fig. 3c) in FREE3 illustrates this problem. The lack of realism of the unresolved mesoscale is even more problematic in WC3. Nonetheless, as shown by Fig. 3 and confirmed by statistical analyses presented in Section 5, FREE3 provides a circulation that captures quite well the mean dynamics of the CCS, which is a critical prerequisite for the application of data assimilation. WC3 is computationally far less expensive than WC10, and allows us to test the efficacy of various data assimilation system configurations, while WC10 is used to evaluate improvements in the circulation dynamics resulting from the best data assimilation configuration.

74

G. Broquet et al. / Dynamics of Atmospheres and Oceans 48 (2009) 69–92

3. The ROMS-IS4DVAR system 3.1. The IS4DVAR algorithm The 4DVAR method was developed by Le Dimet and Talagrand (1986). The ROMS-IS4DVAR system is based on the incremental formulation described by Weaver et al. (2003). We present here the overview and important concepts of this system to aid discussion in later sections, and specific details about the implementation in ROMS can be found in Powell et al. (2008). Let us denote the model state vector as x, composed of the prognostic variables: the sea surface height  at all model grid horizontal locations, and the velocities (u, v), the temperature T and the salinity S at all model grid 3D locations. We represent the NLROMS symbolically as (∂x/∂t) = M(x) where M denotes the model operators. Our aim is to create a circulation estimate {xt }, over a given period [t0 , tf ], that minimizes J NL , the least-squares difference between the available time-sequence of observations {yti } , and {xbt }, i ∈ [1:Nobs ]

the solution from the NLROMS simulation, which is called the background. In the current implementation of ROMS-IS4DVAR, we adjust only the initial-state vector by adding an increment vector ıxt0 to xbt . The boundary conditions, surface forcing, and model parameters are not adjusted. 0

In the incremental approach, it is assumed that {xt } is always close to {xbt }, in which case we seek to minimize the quadratic cost function J = Jb + Jo given by 1 (1) Jb (ıxt0 ) = ıxTt0 B−1 ıxt0 2 Nobs

Jo (ıxt0 ) =

T 1 (Hi Mbt0 →ti ıxt0 − di ) R−1 (Hi Mbt0 →ti ıxt0 − di ) i 2

(2)

1

where Hi is the linearization of Hi , an operator that maps {xti } values onto observation space at ti , Mbt0 →t

is the integration of the tangent linear model (TLROMS) over [t0 , t], di = yti − Hi (xbt ) are the innovai tion vectors, and B and Ri are respectively the model background and observational error covariance matrices. The discussion of how these errors are selected is postponed until Section 4. The integration T

backward in time of the adjoint of TLROMS (ADROMS), Mb t→t0 , yields ∇ J(ıxt0 ), and we seek a solution, ıxat , for which this gradient vanishes and J has its minimum value. A sequence of successive mini0

mizations of J, which is quadratic in ıxt0 , yields a local minimization of J NL = Jb + JoNL , which is difficult

Nobs

to minimize directly, where JoNL = (1/2)

1

(Hi Mt0 →ti (xt0 ) − yti )T R−1 (Hi Mt0 →ti (xt0 ) − yti ). i

Because of the large dimension of the problem (the size of x is O(105 ) in WC3 and O(106 ) in WC10), is identified iteratively by solving a sequence of linear least squares minimizations of (1)+ (2)(so

ıxat 0

called “inner-loops”) repeated with periodic updates of Mb (via so called “outer-loops”, so that J better shapes J NL ). The minimization is achieved using a conjugate-gradient (CG) algorithm based on Fisher (1998). The minimization procedure is terminated when ||∇ J|| ≤ . This typically requires a large number of iterations which is generally not affordable, so instead the number of outer-loops and the number of inner-loops used during each outer-loop were fixed a priori to yield a good estimate of the minimum of J. During each inner-loop, TLROMS is used to propagate the increments forward in time to evaluate (2), and ADROMS yields ∇ J(ıxt0 ) which is used by the CG algorithm to identify the minimum of J. The best estimate for ıxat at the end of the minimization procedure (after Nouter 0 outer . outer-loops) will be denoted ıxN t0 ROMS-IS4DVAR relies on the property of the adjoint to estimate the sensitivities of functions such as J, to changes in the initial condition. If ADROMS indicates that particular features of the circulation are sensitive to changes in the initial conditions, they should be controlled well by IS4DVAR if observed. The theoretical solution minimizing J, as well as the solution provided by the preconditioned CG algorithm, is a correction in the observation space projected into the initial state space by application of the adjoint of the observation operator and the adjoint model, and then smoothed and extrapolated by the covariance B. Following Derber and Bouttier (1999) and Weaver and Courtier (2001) the background error covariance is expressed as B = C, where  is the diagonal matrix of standard deviations, and C is a matrix

G. Broquet et al. / Dynamics of Atmospheres and Oceans 48 (2009) 69–92

75

Fig. 4. A schematic illustration of the assimilation cycle.

of univariate correlations that are modeled as solutions of a diffusion equation. Spatially invariant horizontal and vertical length-scales are used throughout, leading to correlations that are homogeneous and isotropic in both the horizontal and the vertical. While B is univariate, the data assimilation system is effectively multivariate, because the dynamics of ADROMS extrapolate observational information over the whole state-space. As shown in Section 5, this multivariate nature of extrapolations during data assimilation can be useful for generating physically consistent corrections, although the linearity of ADROMS dynamics may generate undesirable features in these extrapolations. 3.2. Assimilation cycle Data assimilation is typically performed over many cycles of length t cycle that span the desired analysis period, thus allowing assimilation windows that are short enough for TLROMS to reliably describe the evolution of increments in NLROMS (see Section 4.1). The procedure consists of the following steps, as illustrated in Fig. 4: • At the beginning of the jth cycle, the initial condition from the previous assimilation cycle, denoted ini,j xt , is used to compute a forecast (hindcast) of length t cycle with NLROMS. This forecast provides j

b,j

the background trajectory {xt }[t ,t

j+1 ]

j

for the jth cycle, as well as an estimate of the optimal model

trajectory when considering only observations available through time to tj . • ROMS-IS4DVAR is applied on the jth cycle [tj , tj+1 ] using Nouter outer-loops of the CG algorithm, to yield an estimate of the optimal increment ıxtNouter ,j for the initial condition at tj . The sum j

Nouter ,j yields xNouter ,j , an estimate of the optimal analysis initial condition. The latter is used xini t + ıxt t j

j

j

to compute the analysis {xtNouter ,j }[t ,t j

j+1 ]

with NLROMS during the jth cycle, which represents the esti-

mate of the optimal model trajectory when considering observations available until tj+1 . The state ini,j+1

xtNouter ,j is then used as the starting point for the (j + 1)th cycle: xt j+1

j+1

= xtNouter ,j . j+1

4. Set up of the data assimilation experiments The IS4DVAR set-up was chosen to be consistent between WC3 and WC10, so that conclusions based on sensitivity studies in WC3 provide guidance for the computationally more expensive WC10.

76

G. Broquet et al. / Dynamics of Atmospheres and Oceans 48 (2009) 69–92

Discussions of specific sensitivity studies will be postponed until Sections 5.2 and 5.3, but are referred to in subsequent sections. 4.1. ROMS-IS4DVAR parameters The chosen configuration of ROMS-IS4DVAR addresses several critical points: (i) The length of each assimilation cycle must ensure that the tangent linear assumption (the assumption that the time evolution of the increments obeys the tangent linear model dynamics) that yields the linearized problem in (1) and (2) is approximately valid. However, the assimilation window must be long enough to span periods with sufficient observations capable of constraining the model dynamics. Veneziani et al. (2009b) have shown that in WC10, for perturbations typical of those encountered in IS4DVAR, the linear assumption is approximately valid for ∼ 14 days, which is our choice for the data assimilation window (this assumption being a priori also valid at the lower resolution of WC3). Additional experiments (not presented here) conducted in WC3 with a 7 days assimilation cycle revealed that the assimilation and forecast results are not sensitive to the assimilation window in this time range. (ii) The number of inner- and outer-loops used is a trade-off between computational cost, the level of convergence of the CG algorithm, and the need to update the cost function shape. Powell et al. (2008) used 14 days assimilation windows in the IAS and found that using several outer-loops had a weak influence on the cost function convergence, in that for a given total number of loops (i.e. the product of the number of inner-loops and outer-loops), the final cost function was similar irrespective of the combination of inner- and outer-loops. Powell et al. (2008) also found that ∼ 20 loops in total (inner- times outer-loops) was sufficient to yield an acceptable estimate for the minimum of J, in the sense that the rate of change in J indicated that differences between its final value and the optimal solution are small compared with the difference between the initial and final values. Similar results were found using WC3, which was finally run with 1 outer-loop and 21 inner-loops, although, the computational cost of WC10 limited the total number of loops to 11. The issue of the rate of convergence to the minimum of J is studied further in Section 5. (iii) The parameterization of B significantly influences the way that observational information is extrapolated to the non-observed variables, and smoothed over the observed variables. However background errors are very poorly known. The standard deviations  dictate the relative weights of confidence between the different components of the background initial state in Jb and the observations in Jo . Monthly estimates of  for WC3 (WC10) were calculated with FREE3 (FREE10) during 1999–2004. This should yield a reasonable representation of , since Veneziani et al. (2009a) showed that model variability in EKE, SST and hydrography are in good agreement with observations. The main source of background error is linked to model bias and to position errors of circulation features in the model (which decrease the correlations between the model and the observations). The dynamical estimate of  combines with the dynamical spatial extrapolation of the adjoint to produce a spatial extrapolation of observational information that is related to the coastal ocean dynamics. Sensitivity tests were performed by rescaling  in order to explore the effects of modifying relative contributions of B and Ri to J. However, the best results have been obtained when using  derived directly from model intrinsic variability. The decorrelation length scales of C should be similar to the natural decorrelation length scales of the related variables, and should be linked to the control of important dynamical features. Increasing the decorrelation lengths allows for greater extrapolation of the information from sparse data, but can overly smooth the information from dense observing systems. Identical decorrelation lengths were used in WC3 and WC10. Preliminary experiments showed that with the nearly 70 km spacing between stations and lines in CalCOFI data (see below), under the constraint of using homogeneous values, it seems appropriate to use a 50 km horizontal decorrelation length scale for all variables (we impose the same error space scales on all variables). This is more than the spacing of SSH and SST available data but is close to the resolution of WC3, which means the horizontal extrapolation of data is weak. In the vertical, a 30 m decorrelation scale prevents surface

G. Broquet et al. / Dynamics of Atmospheres and Oceans 48 (2009) 69–92

77

data from directly exerting an influence well below the surface. Experiments with greater horizontal and vertical decorrelation length scales in WC3 indicate that solutions are overly smoothed when using longer correlation scales. (iv) Estimates of Ri are based on a diagonal parameterization, where the errors of different observations are assumed to be uncorrelated. This may be a poor assumption for gridded products such as used here, but accounting for correlated observational errors raises some theoretical and technical difficulties. The variances of these errors should reflect both the measurement errors and errors of representativeness. The latter are not well known, although variations in this term can be used as an alternative way of controlling the contribution of the model background and the observations to J. It also allows the various types of observations to be weighted differently. For example, imposing more accuracy on in situ data that are very few in number compared to satellite data can be useful. Preliminary experiments in both WC3 and WC10 have shown that decreasing the observational error assigned to in situ data leads to significant improvements at depth, so relatively low measurement errors are used here (see below Section 4.2). The variance of the measurement errors, o , depends on the physical nature of the observations and their origin. The standard deviation, g , of the observations in a grid cell is used as an indicator for the error of representativeness. For this experiment, the observational error is assigned to either the measurement error or the error of representativeness, whichever is higher. 4.2. Observations for data assimilation and data/model comparisons Most observations available for the CCS during the period 1999–2004 are satellite remotely sensed surface data. The surface datasets used here are • Aviso SSH anomalies: the product used is actually a merging of Ssalto-Duacs data (from TOPEX/Poseidon, Jason-1, Envisat and GFO measurements), provided by Aviso and the Centre National d’Etudes Spatiales (CNES). The gridded data are useful for resolving the dynamically important SSH gradients. The anomalies, relative to the CNES dynamic SSH mean, are converted into total dynamic height consistent with the model, using the process described in Powell et al. (2008). The steric signal is removed (since it is not resolved by ROMS) using the Willis et al. (2004) database, and a mean dynamic height of the model over 1999–2004 is used to compute the total dynamic height. The data are available every week with approximately a 1/3◦ resolution in the CCS. An observation error standard deviation o = 2 cm is used for SSH measurements, as it is the usual measurement precision associated with altimetric data. • SST from COAMPS: this product is based on an OI system, Navy Coupled Ocean Data Assimilation (NCODA; Cummings, 2005), that can be executed in two-dimensional mode to provide SST and sea ice concentration lower-boundary conditions. The ocean data types assimilated for the 2D analysis include: remotely sensed SST, sea-ice concentration and in situ surface observations from ships and buoys. An ocean data quality-control system is fully integrated as well. This provides SST data set each day, at a spatial resolution higher than that of WC3 and WC10 in the CCS region, and which is consistent with the COAMPS winds and other surface fields. In order to use no more than one datum per horizontal grid point in each configuration, the data were interpolated on each model grid (these are the only data that were interpolated on the model grid). For the SST measurements o = 0.4 ◦ C. In situ data consist of hydrographic data (T and S profiles) of CTD measurements from: • The CalCOFI cruise lines: the T/S data used were collected along six cruise lines conducted four times every year in the SCB. The measurements extend to 500 m depth, and the stations are ∼70 km apart. • The GLOBal ocean ECosystems dynamics (GLOBEC), Northeast Pacific (NEP)/CCS mesoscale survey cruises (Batchelder et al., 2002). Up to five cruise lines are conducted four or five times every year

78

G. Broquet et al. / Dynamics of Atmospheres and Oceans 48 (2009) 69–92

Fig. 5. (*) CalCOFI, (◦) GLOBEC and (+) NMFS typical locations of measurement used.

off the coast of Oregon, with observations of T/S extending to more than 1000 m depth. The space between measurement stations ranges from 5 km to 30 km. • The NOAA National Marine Fisheries Service (NMFS) datasets (Baltz et al., 2006): consisting of closely spaced measurement stations in the vicinity of Monterey Bay and available from May to June every year, providing T/S observations vertically every 2 m to greater than 500 m depth. The typical sampling arrays for these measurements are shown in Fig. 5. We use o = 0.1 ◦ C (0.01 psu) for T (S) data which is probably optimistic, but serves to increase the weight of in situ data relative to the satellite data, especially SST, during the assimilation process. 4.3. Experimental protocol Simulations on WC3 and WC10 were initialized respectively from FREE3 and FREE10 on January 1st 1999 or 2000, and data assimilation performed over the period 1999–2004 or 2000–2004. Several types of data assimilation experiments were conducted, using different subsets of the observational data described above. 5. General results from surface and in situ data assimilation The results presented here explore the effects of simultaneously assimilating SSH, SST and vertical T/Sprofiles. They also demonstrate how the data can combine to improve the model circulation estimates, and raise some interesting issues about the optimality of this combination. It should be noted that WC3 (WC10) uses the k–  (k– ω) GLS vertical mixing parameterization. The impact of vertical mixing parameterizations on the assimilation results will be considered in Section 5.3. Variability in SST, SSH and in T/S profiles along the coast result from local and remote forcing mechanisms, and using data assimilation we seek to improve these aspects of the circulation. Veneziani

G. Broquet et al. / Dynamics of Atmospheres and Oceans 48 (2009) 69–92

79

Table 1 Main experiments conducted in the CCS. Reference

Configuration

Data assimilated

Period

FREE3 FREE10 ALL3 ALL10 CTD3 SAT3

WC3 WC10 WC3 WC10 WC3 WC3

No data assimilation No data assimilation SSH, SST, CalCOFI and GLOBEC data SSH, SST, CalCOFI and GLOBEC data CalCOFI and GLOBEC data SSH and SST data

1999–2004 1999–2004 2000–2004 1999–2004 2000–2004 2000–2004

et al. (2009b) showed that over a 14-day period, the circulation in the vicinity of Monterey Bay is about equally sensitive to variations in the model three-dimensional initial conditions and surface forcing. Therefore, over short assimilation periods, even in the absence of corrections to the surface forcing, the initial conditions are expected to exert an important control on the coastal dynamics of the CCS. These sensitivities indicate that there will be a seasonal dependence of the skill of ROMS-IS4DVAR. Statistics of the differences between the model and the observations are studied as quantitative measures of the improvement in the model circulation. These statistics are estimated spatially and temporally using various observational data sets {yj }j ∈ [1...L] and model solutions {xj }j ∈ [1...L] interpolated to the observation times and locations. They describe the error through the separation of the Root Mean

L

Square (RMS) error ((1/L)

j=1

L

(xj − yj )2 )

1/2

into bias given by (1/L) 2 1/2

¯ ) ation (STD) given by ((1/L) j=1 (xj − x¯ − yj + y) and according to (Taylor, 2001):

L

where x¯ = (1/L)

j=1

L

(RMS Error)2 = (Error BIAS)2 + (Error STD)2

(xj − yj ) and standard devi-

j=1

(xj ) and y¯ = (1/L)

L

j=1

(yj ),

(3)

IS4DVAR is formulated based on the assumption that the model, observations and resulting analyses are unbiased. Biases in the model dynamics and forcing are unavoidable but are interpreted by the assimilation system as part of the error STD which renders the resulting analysis suboptimal. The error STD can also be written in terms of the correlations between model and the data, the L ¯ where referred to here as the model data correlation (MDC) given by (1/Lx y ) j=1 (xj − x¯ )(yj − y)

L

1/2

L

1/2

¯ 2 ) , and the standard deviation error (SDE) and y = ((1/L) j=1 (yj − y) x = ((1/L) j=1 (xj − x¯ )2 ) which is x − y , the difference between the Model STD and the Obs STD: (Error STD)2 = (SDE)2 + 2 × Model STD × Obs STD × (1 − MDC)

(4)

The SDE is indicative of the level of agreement between circulation features in the model and in data, while a high MDC indicates that, when the SDE is low, these features occur at similar locations and times in both the model and the observations. 5.1. Assimilation of surface and in situ data We first consider the results from the assimilation of all available data presented in Section 4 except NMFS data. NMFS data were excluded for independent evaluations of the resulting analyses. Improvements in the circulation that are quantified using non assimilated data confirm a positive impact of data assimilation on the model circulation as a whole. The assimilation of both surface and in situ data provides the best estimate of the circulation in the CCS. However, there is some competition during the assimilation process between surface and in situ data which will be analyzed via separate experiments in Section 5.2. This competition is related to limitations of the model and data assimilation system. Experiments using SSH, SST, CalCOFI and GLOBEC data assimilation in WC3 and WC10 are referred to respectively as ALL3 (for the period 2000–2004) and ALL10 (for 1999–2004). The notation used for each experiment is summarized in Table 1.

80

G. Broquet et al. / Dynamics of Atmospheres and Oceans 48 (2009) 69–92

Fig. 6. Evolution of (a and c) ||∇ J||/||∇ J||inner-loop 1 and (b and d) J (continuous line) and J NL (crosses) during the assimilation cycle t /2. from March 25th to April 8th 2000 in (a and b) ALL3 and (c and d) ALL10. Horizontal thick dotted line in (b and d) is Nobs

5.1.1. Cost function behavior The decrease of J towards a minimum value during a typical assimilation cycle for ALL3 and ALL10, is illustrated in Fig. 6 for the period 25 March–8 April 2000. J along with J NL (Fig. 6b and d) and ||∇ J||/||∇ J||inner-loop 1 (Fig. 6a and c) (estimated using respectively TLROMS and ADROMS) are shown as functions of the number inner-loops. The largest change in J and decrease in the gradient norm occur during the first five inner-loops. The relative decrease in the gradient norm from one iteration changes from ∼50% between the two first inner loops to less than 20% after the third inner loop. This indicates that 10 inner-loops are sufficient to reach an acceptable solution in both ALL3 and ALL10. During a given cycle of assimilation, if B and Ri were the “true” Gaussian error covariances of the t /2, with background and observations, the theoretical minimum of J would have a mean Jmin = Nobs



t /2 (Weaver et al., 2003), where N t standard deviation of Nobs is the total number of assimilated obs t ∼336, 000 for ALL10 per cycle. Fig. 6b and d indit observations. Here Nobs ∼34, 000 for ALL3 and Nobs cates that during the assimilation cycle 25 March–8 April 2000, J converges to ∼1.2Jmin in ALL10 and ∼1.4Jmin in ALL3, further from the theoretical value in the latter case. This illustrates the result that the means of J/Jmin at the last inner loop, over every assimilation cycles in ALL10 and ALL3, are respec t t /2∼180 for ALL3 and Nobs Nobs /2∼580 for ALL10 per cycle so the final tively ∼1.25 and ∼1.35. values for J are relatively far from Jmin when compared with the standard deviations of the theoretical distributions for the minimum of J. However, the differences between final values for J and Jmin are relatively small when compared to the differences between initial values for J and Jmin . They are also relatively small considering that the theoretical distribution for the minimum of J is based on a perfect

G. Broquet et al. / Dynamics of Atmospheres and Oceans 48 (2009) 69–92

81

parameterization of error covariances, for which uncertainties may lead to changes in J far greater than the differences (J − Jmin ) noted above. Even though the error covariances are a little underestimated for WC3 (the error of representativeness is likely an underestimate at low resolution), our choice of the relative weighting in observational and background errors appears justified, especially for WC10. However, the analysis of Jo and Jb (not shown) indicate that during ALL3 (ALL10) Jo ∼10Jb (∼103 Jb ) and dominates J, suggesting that the background error covariances seem to be overestimated while the observational error covariances seem underestimated. As a result Jb continues to increase even after J has reached a value that appears close to its minimum, which implies that xt can deviate from xbt without changing the overall fit of the model to the observations. These adjustments were found to occur at scales that decrease rapidly during the first few inner loops, and after 20 inner loops, the increments are characterized by very localised small scales which appear spurious and undesirable. To remedy the growth of Jb we could increase Ri , but preliminary experiments have shown that this prevents the model from sufficiently fitting the observations. Alternatively, terminating the minimization algorithm after 10 or 20 inner-loops prevents the model from diverging too much while ensuring a satisfactory fit to observations. Lowering  was tested, leading to better convergence of Jb but also to a significant increase in the asymptotic value of Jo . The analyses of Jo also reveal the reduction of the contributions of in situ S and T during the inner-loops, indicating that in situ data were weighted sufficiently to impact J despite the fact that there are far fewer in situ observations than satellite data. This is a second justification of our choice of observational errors. J NL (Fig. 6b and d) can only be computed in the outer-loops, so only two values are reported. The final value of J NL is larger than the final value of J, which is typical in 4DVAR (Tshimanga et al., 2008), but regardless J NL is still reduced considerably by data assimilation. 5.1.2. Error statistics at surface The improvements in the model circulation associated with data assimilation are illustrated in Fig. 7, which shows the evolution of each term in (3) and (4) for SSH and SST, during each assimilation cycle for ALL10 (considering all observations available during each 14 days cycle). Fig. 7 shows the error of 14 days forecasts initialized from the end of each analysis cycle. Results during the forecast period yield an estimate of the improvements arising from past analyses. For FREE10, the RMS error for SSH increases slightly over time (Fig. 7a). The RMS error is composed of a bias (Fig. 7b) that varies seasonally, becoming more negative over time, and an error STD (Fig. 7c) with similar amplitude and a slightly increasing trend. The latter is mainly characterized by a small MDC (Fig. 7e) since the SDE is small (Fig. 7d), indicating that the model produces realistic SSH structures, but position errors are present. For SST, the RMS error for FREE10 (Fig. 7f), also increasing slightly, is composed of a negative bias (Fig. 7g) and of an error STD with comparable amplitude (Fig. 7h). Both the bias and error STD vary on seasonal and interannual timescales. The SDE of SST (Fig. 7i) is large, and mainly negative (positive) during fall–winter (spring–summer), showing that the model tends to enhance the natural seasonal variability. In other words, the SSH and SST Model STD increase too much during upwelling events and the periods of highest eddy activity, and decrease too much during the winter. The bias on the SSH also intensifies the natural SSH tendencies near the coast, with the model overestimating (underestimating) SSH during winter (summer). FREE3 exhibits similar behavior with comparable errors. Fig. 7b and g shows that even though the bias in SSH and SST is not formally accounted for in the formulation of IS4DVAR, data assimilation is very effective at reducing them. In fact, the system displays an ability to control all characteristics of the error, as evidenced by the reduction in SDE (Fig. 7d and i) and the increase in MDC (Fig. 7e and j). Forecasts from ALL10 sometimes exhibit errors as large as those of FREE10 but in general forecast errors are smaller, indicating that the adjustments made during each analysis lead to improvements over several cycles, despite the fact that a portion of the error is associated with errors in the surface forcing, which is not adjusted. During the assimilation, RMS error is generally reduced by more than a factor of 2 in SST (Fig. 7f), at both the analysis and forecast times. However, the skill of the data assimilation system varies seasonally and experiences an increase in error during spring and summer, probably because the sensitivity of the model to atmospheric forcing is maximal during this period (Moore et al., 2009; Veneziani et al., 2009b). For SSH, the RMS error in ALL10 is reduced less when compared to that of FREE10 (Fig. 7a)

82

G. Broquet et al. / Dynamics of Atmospheres and Oceans 48 (2009) 69–92

Fig. 7. Error statistics from Eqs. (3) and (4) for each assimilation cycle in WC10 for (a–e) SSH and (f–j) SST, relative to all Aviso and COAMPS data available during each cycle. Light grey: FREE10; dark grey: forecast in ALL10; black: analysis in ALL10.

than for SST, especially because the MDC has not improved much in ALL10 (Fig. 7e). This can be related to the relatively small number of SSH observations when compared to the number of SST observations available for WC10 (the ratio being smaller in WC10 than in WC3) and to the difficulties in recovering SSH gradients from their relatively sparse distribution at such a resolution, which explains that in ALL3, errors in SSH are clearly reduced more when compared to FREE3 (not shown). Another reason for the lower decrease of the error in SSH is that a significant fraction of the SSH information is lost in the model via geostrophic adjustment in the form of barotropic waves as a result of the imbalances that exist in the assimilation increments between SSH, T and S. Based on the balances imposed by the adjoint model, the main adjustments to fit SSH appear in the T and S initial fields. Changes in the T and S initial fields are also related to the assimilation of SST, and in situ T and S data that do not lead to necessarily compatible corrections, due to the limitation of the system, and to the fact that the mean dynamic height of the SSH observations is not prescribed with accuracy. 5.1.3. Horizontal circulation structure The mean June SST and summer mean SSH for 2003–2004 for ALL3 and ALL10 are shown in Fig. 8. Comparing Fig. 8 with Fig. 3 illustrates the corrections made by IS4DVAR to SSH and SST. Key features to note are changes in the position of the CC and of the coastal jet and a quantitatively altered meandering of the coastal jet. These are key points considering that no adjustments to the forcing were made. Offshore extensions of SST features are suppressed by assimilation, especially in WC10 (compare Fig. 3a, b with Fig. 8a and b). The erroneous dome of SSH near the NW open boundary in WC3 (Fig. 3c) is eliminated by IS4DVAR (Fig. 8c). Even if the open boundary conditions are not adjusted, the data assimilation yields fields inside the domain of both WC3 and WC10 that are far more coherent with the ECCO-GODAE boundary conditions, and thus the model fields have reduced boundary artifacts. However, some fields do possess discontinuities that, in ALL3 as in FREE3, generate an anticyclonic

G. Broquet et al. / Dynamics of Atmospheres and Oceans 48 (2009) 69–92

83

Fig. 8. Same as Fig. 3 but for (a and c) ALL3 and (b and d) ALL10.

structure along the western and northern boundaries. During winter (not shown), the mean SSH is reduced. 5.1.4. Vertical structure of the circulation The subsurface ocean structure is also strongly affected by data assimilation. Despite an horizontally sparse distribution of in situ data, data assimilation yields smooth fields on the analysis initial conditions (not shown) because of the dynamical influence of the adjoint model spread and smoothing by B (even if the latter is applied with small length scales). Fig. 9 shows subsurface error statistics for FREE10, FREE3, ALL10 and ALL3 using CalCOFI (Fig. 9a–h) and NMFS measurements (Fig. 9i–l) during period 2000–2004, and illustrates the typical behavior of WC10 or WC3. For FREE10 (and FREE3) in winter, the bias in T is positive between 50 and 150 m depth, and maximum around 100 m depth, although slightly negative near the surface (Fig. 9a). During summer, the bias is positive above 200 m and decreases with depth (Fig. 9e). The error STD is also maximum around 100 m (50 m) depth during winter (summer) (Fig. 9a and e). The Model STD, Obs STD, and the structure of the bias suggest that in winter the error is maximum at the model thermocline (Fig. 9a), which is deeper than observed. In summer the model thermocline is comparable in depth to the observations, but has a lower STD (Fig. 9e). For subsurface S, there is a negative bias year-round centered at 250 m depth and an error STD maximum in the upper 100 m depth of the water column (Fig. 9b and f). The statistics for ALL10 and ALL3 in Fig. 9c, d and g, h show that the vertical structure has been corrected well, especially in S. Statistics computed with respect to GLOBEC data are similar to those

84

G. Broquet et al. / Dynamics of Atmospheres and Oceans 48 (2009) 69–92

Fig. 9. Statistics of subsurface errors in (a, c, e, g, i and k) T and (b, d, f, h, j and l) S computed from (a–d) CalCOFI measurements in January–February, (e–h) CalCOFI measurements in June–July and (i–l) NMFS measurements in May–June during the period 2000–2004 as a function of depth for (a, b, i and j) FREE10 (e and f) FREE3 (c, d, k and l) ALL10 and (g and h) ALL3: Error Bias (black), Error STD (dark grey), Model STD (grey) and Obs STD (light grey).

with respect to CalCOFI (not shown). Subsurface warm biases in winter are larger in ALL10 when compared to FREE10 (Fig. 9a and c), a behavior that is linked to the assimilation of surface data and that also limits improvements in summer and in ALL3 (see discussions on this point in Section 5.2). However, warm biases during summer in WC10 and WC3 and warm biases during winter in WC3 become weak cold biases after data assimilation (Fig. 9e and g). Very good agreement between the

G. Broquet et al. / Dynamics of Atmospheres and Oceans 48 (2009) 69–92

85

Model STD and Obs STD for T is reached (Fig. 9c and g), and the model reproduces the observed depth of thermocline. The error STD remains similar below 200 m, maximum around 100 m in winter and decreases with increasing depth during summer, but a net correction of nearly 0.3 ◦ C is reached in the upper 150 m (Fig. 9c and g). Error STD values in ALL10 or ALL3 are generally less than 1.2 ◦ C. For S, the negative bias that occurs throughout the year is replaced by a positive bias near the thermocline, and is completely removed elsewhere (Fig. 9b, d and f, h). The very good correction of the SDE (shown by the better fit of the Model STD to the Obs STD) and the clear improvement of the error STD may be due to the very low observational error associated with S data, but based on these results, the model, even with coarse resolution, is able to represent very well the water column structure. Fig. 9i–l shows independent error statistics that are based on the unassimilated NMFS data for ALL10 which are similar for ALL3 (not shown). Fig. 9i–l indicates that the assimilation of CalCOFI and GLOBEC data has an important and positive impact on the mean T and S structure far from the observation locations. Note that this really is a remotely influenced effect and not a local SST and SSH effect (as demonstrated by experiments with the assimilation of in situ data or surface data alone, not shown). This is a key result because data assimilation systems often experience difficulties in extrapolating information from assimilated data. The dynamical part of this extrapolation, achieved by ADROMS, leads to a general adjustment of the mean water column near the coast. The role played by the spatial extrapolation of data by B is secondary because short horizontal length scales were used to model C. Dynamical propagation through NLROMS of IS4DVAR increments from the CalCOFI and GLOBEC cruise tracks is also likely to play a role in the correction of the circulation in the vicinity of Monterey Bay. Data assimilation significantly reduces the bias relative to the NMFS salinity (Fig. 9j and l), and yields a better fit of the Model STD to the Obs STD for S. Despite a correction of the cold bias below 150 m, data assimilation does not improve T in the upper 150 m (Fig. 9i and k). In fact, assimilation of surface data has a somewhat negative impact on the bias in T (see Section 5.2). However WC3 experiments with the assimilation of CalCOFI and GLOBEC data alone leads to a reduction of the bias relative to the NMFS temperature (not shown). The adjustments in SSH and on the structure of T and S, even far from the in situ data locations, and the fact that the model dynamics adjust the velocity fields via ADROMS and geostrophic adjustment, lead to significant modifications of the currents in the whole CCS, and to adjustments of the upper ocean baroclinic structure, especially in coastal regions. To illustrate, Fig. 10 shows the summer mean meridional velocity along section at 37◦ N during 2003–2004 from ALL3 and ALL10. Comparing Fig. 10 with Fig. 2 reveals that data assimilation has modified the CC and the CUC, and the poleward current abutting the CUC is weakened considerably.

5.2. Assimilation of in situ data and surface data separately Exploring the effects of assimilating SSH and SST or vertical T/S profiles separately demonstrates the relative contributions of these data during ALL3 and ALL10. However, as detailed in this section, it also demonstrates that improvements from assimilation of in situ data can be degraded by the assimilation of surface data. In WC3, experiments using SSH and SST or CalCOFI and GLOBEC data assimilation separately are referred to as SAT3 and CTD3 respectively (cf. Table 1). Both were performed during the period 2000–2004. Comparable experiments in WC10 were only performed during 2000, and exhibited very similar behavior to WC3, so only the results for CTD3 and SAT3 are presented. In situ data assimilation has a quantifiable influence on surface fields, especially on SSH near the coasts. Fig. 11, which should be compared with Fig. 3c, d, and Fig. 8c, d, indicates a relocation of the coastal features in FREE3 (Fig. 3c) to CTD3 (Fig. 11a). CTD3 (Fig. 11a) reproduces some features of SAT3 (Fig. 11b) when compared to FREE3 (Fig. 3c), suppressing the spurious circulation feature in the north–west corner of WC3. However, in CTD3 the coastal jet is more deviated from the coast south Cape Mendocino than in FREE3 while in SAT3 it gets closer to the coast than in FREE3. The similarity of surface fields and error statistics for SSH and SST in ALL3 (Fig. 8c) and SAT3 (Fig. 11b) shows that the influence of SSH and SST data assimilation dominates the influence of in situ data. Along the CalCOFI and GLOBEC cruise tracks, the circulation patterns in ALL3 are in better agreement with SAT3 than with CTD3.

86

G. Broquet et al. / Dynamics of Atmospheres and Oceans 48 (2009) 69–92

Fig. 10. Same as Fig. 2 but for (a) ALL3 and (b) ALL10.

Fig. 12 shows error statistics computed from T and S CalCOFI data for CTD3 and SAT3 and should be compared with Fig. 9e–h. The T and S error STD is generally reduced by assimilation of SSH and SST especially near the surface. Surface data assimilation also increases the mean salinity at all depths which generates a strong salinity bias in the upper 100 m (that does not persist when assimilating all data), but decreases the negative bias below 150 (Fig. 9f, h and d). However, the largest impact of surface data assimilation is on the bias in T which increases during winter and summer as evidenced by a more diffuse thermocline (Fig. 9e, g and c). Near the surface, the cold winter bias is replaced by a warm bias during the assimilation of SST leading to an overall heating of the upper water column (not shown for SAT3 but see Fig. 9a and c). On the other hand, assimilation of in situ data (Fig. 12a and b) removes the warm bias of winter and summer replacing them with a weak cold bias with errors on T (S) ≤1.5 ◦ C (0.3 psu), regardless of depth. The increase of the bias in T near the thermocline relative to

Fig. 11. Summer SSH mean in m for the period 2000–2004 in (a) CTD3 and (b) SAT3.

G. Broquet et al. / Dynamics of Atmospheres and Oceans 48 (2009) 69–92

87

Fig. 12. Statistics of subsurface errors in (a and c) T and (b and d) S computed from CalCOFI measurements in June–July during the period 2000–2004 as a function of depth for (a and b) CTD3 and (c and d) SAT3: Error Bias (black), Error STD (dark grey), Model STD (grey) and Obs STD (light grey).

the winter CalCOFI data (Fig. 9a and c) or the spring NMFS data (Fig. 9i and k) for ALL10 (discussed in Section 5.1.4) is therefore associated with the assimilation of SSH and SST data. Fig. 13 shows comparisons of observed T along the CalCOFI line 77 in July 2002 with similar sections from FREE3, SAT3 and CTD3, and illustrates the weakening of the model thermocline associated with the assimilation of SST and SSH (Fig. 13b and c). Fig. 13d shows how the adjustments made during CTD3 reproduce the observed vertical structure (Fig. 13a). Not only is the thermocline in CTD3 less diffuse, but the isotherms slopes near the coast are more realistic and indicative of the signatures on T and S of poleward subsurface currents. Weakening of the thermocline and increase of the bias in T during the assimilation of SST and SSH (Fig. 13c) is related to the extrapolation of surface observations via the adjoint dynamics and the vertical diffusion operator of C. This could have been expected because the vertical diffusion operator introduces vertical corrections that are too large in the thermocline. The sensitivity of this warm bias to the vertical extrapolation of surface observations is studied in Section 5.3. The best results when comparing the model and in situ data sets (the assimilated CalCOFI and GLOBEC data or the non assimilated NMFS data) are obtained when assimilating only in situ data. Fig. 14 shows the RMS error relative to all CalCOFI and GLOBEC data available per 14 days cycle, during FREE3 and CTD3. Typically, the RMS error in T and S is reduced by more than a factor of 2 by data assimilation, confirming the good agreement between the observations and the analysis. However, the forecast error lies closer to that of the free simulation than to the analysis, indicating that the forecast generally degrades quickly. The loss of accuracy at depth when comparing CTD3 and ALL3 or when comparing analysis and forecast in CTD3 may be linked to the configuration of ROMS-IS4DVAR, such as the horizontal and vertical correlation lengths implicit in C. Its potential link to the duration of the forecast was explored by using 7 days assimilation cycle, without any improvement (see Section 4.1i). 5.3. Background error vertical correlations and model vertical mixing The negative influence of surface data on subsurface T discussed in Sections 5.1 and 5.2 is linked to the vertical extrapolation of SST and SSH observations by IS4DVAR, and depends primarily on: (i) the diffusion operators associated with the vertical correlations of C, and (ii) the vertical diffusion parameterization in NLROMS. (i) Preliminary experiments revealed that vertical correlations > 30 m are probably too large close to the surface, and impact substantially the subsurface structure. At the same time, deep terrainfollowing coordinates off the shelf are separated by distances >60 m below about 200 m depth.

88

G. Broquet et al. / Dynamics of Atmospheres and Oceans 48 (2009) 69–92

Fig. 13. T along the CalCOFI line 77 in July 2002 (white points indicate CalCOFI stations at standard levels). (a) CalCOFI data; (b) FREE3; (c) SAT3; (d) CTD3.

G. Broquet et al. / Dynamics of Atmospheres and Oceans 48 (2009) 69–92

89

Fig. 14. RMS Error for each assimilation cycle in WC3 for (a) T and (b) S relative to all the CalCOFI and GLOBEC data available during each cycle (the data availability is not continuous in time). Light grey: FREE3; dark grey: forecast in CTD3; black: analysis in CTD3. Horizontal doted lines: global RMS Errors considering all data available during the period 2000–2004.

Thus a 30 m or less decorrelation scale leads to vanishing correlations in this region. As shown by Weaver et al. (2003), allowing the vertical correlations to vary with depth is important. Experiments conducted by setting to zero all vertical correlations in C are studied here. (ii) The experiments presented in Section 5 were conducted in WC3 using the k–  GLS vertical mixing scheme. Because of the issues associated with the vertical extrapolation of the observations, experiments were also performed using the k– ω GLS and the KPP vertical mixing schemes. Fig. 15 shows the bias in T, relative to the CalCOFI data in winter during the period 2000–2004, for FREE3, SAT3 and CTD3, with or without vertical correlations in C, for several vertical mixing schemes. The bias is systematically lower when the vertical correlations in C are set to zero. This confirms that the use of diffusion operators for C is problematic. It is less sensitive in CTD3 (Fig. 15c) as the assimilated data have a higher vertical resolution than the model, so that the vertical correlations in C only smooth the corrections but do not extrapolate them to non-observed parts of the water column. Each of the vertical mixing parameterizations for NLROMS have their own advantages and drawbacks. The bias when using k– ω (not shown) or k–  GLS vertical mixing parameterization is very similar. The KPP mixing reduces the warm bias in FREE3 and SAT3 in the first 150 m depth (Fig. 15a and b), but generates a cold bias in the latter below and increases the cold bias in CTD3 (Fig. 15c). Considerable improvements when assimilating surface data are evident in SAT3 when suppressing vertical correlations in C and using the KPP vertical mixing. This confirms that the bias in T is linked to the configuration of ROMS-

Fig. 15. Subsurface error biases in T computed from CalCOFI measurements from January to February during the period 2000–2004 as a function of depth for (a) FREE3, (b) SAT3 and (c) CTD3, using the k- GLS (black) or KPP (grey) vertical mixing parameterization. For SAT3 and CTD3, continuous lines designates the absence of vertical correlations in C and dashed lines the use of vertical correlations in C.

90

G. Broquet et al. / Dynamics of Atmospheres and Oceans 48 (2009) 69–92

IS4DVAR. Reducing the vertical correlations in C limits the extrapolation of data instead of improving it, which has a positive impact because the vertical resolution in the data is higher than the vertical resolution of WC3. Therefore vertical correlations were used in ALL10. Vertical correlations in nature are nonzero, and improvements in the estimation of C are needed.

6. Conclusion This paper presents an application of ROMS-IS4DVAR to the CCS using real observations. This first application of a 4DVAR system to control a fully realistic configuration of the whole CCS was successful in both the 1/10◦ resolution configuration and the coarse 1/3◦ configuration used for preliminary sensitivity tests. Despite competition between the different types of observational data during the assimilation process, the system allows for a coherent assimilation of all available data. Specifically:

(i) The difference between the model and either in situ and surface assimilated data was reduced by more than a factor of 2 compared to the case without data assimilation. The already good agreement of the model simulation without data assimilation to the observations was improved, leading to a more realistic representation of the circulation in the 1/10◦ model. (ii) The difference between the model and independent observations near Monterey Bay that were not assimilated into the model demonstrates a sensible extrapolation by the model dynamics of the information from observed regions to the remote portions of the circulation. (iii) The difference between the model and observations during forecasts, initialized from the end of the analysis in the previous assimilation cycle, indicate that information from each assimilation cycle still leads to improvements in the model circulation for at least 14 days. (iv) The entire circulation of the CCS appears to be sensibly adjusted despite the limited temporal and spatial coverage of the observations.

This has been facilitated by a sensible choice of IS4DVAR configuration, and by processing of the observations prior to assimilation. The role of the dynamical inversion from the use of ADROMS was also shown to be critical, allowing for sensible multivariate adjustments that were necessary to adjust the SSH fields, and for limiting the influence of imperfections in C which is isotropic, homogeneous and univariate. The ability of the adjoint model to extrapolate adjustments that are consistent with dynamics of the model demonstrates the strength of adjoint methods. Further improvements of the CCS ROMS-IS4DVAR system are anticipated that include adjustments to the surface forcing and open boundary conditions in addition to the model initial conditions. The adjustment of initial conditions that may be at odds with open boundary conditions and surface forcing is problematic, and the role of surface forcing, especially the wind, is critical in a region such as the CCS. The initial conditions control only a part of the circulation in the model, and controls by the surface forcing are likely to be important during certain times of the year, such as the summer when a reduction of skill was found when adjusting only the initial conditions. The development of a multivariate B should also increase the ability to extend adjustments to the whole state vector in a sensible manner. Known T–S relations in particular have been shown by Ricci et al. (2005) to improve results at depth. Performance of the system are also shown to be sensitivitive to the choice of . Flowdependant updates of  and B are important as demonstrated in ROMS by Powell et al. (2008) and could strongly improve performance. Finally, better preconditioning of the minimization algorithm, based on estimates of the Hessian of J, are under development. This should modify the problematic behavior of the lack of sensitivity of the cost function to the increments during the minimization after ∼ 10 inner loops. Nonetheless, despite the limitations and issues with the current system, the prospects for an operational forecast and analysis system based on ROMS-IS4DVAR are bright, and will pave the way for further scientific understanding.

G. Broquet et al. / Dynamics of Atmospheres and Oceans 48 (2009) 69–92

91

Acknowledgements This work was supported by the National Oceanographic Partnership Program (NOPP) project NA05NOS4731242. We would like to thank David Foley (NOAA) and Patrick Heimbach (MIT) for useful scientific discussions and for providing data used in this study, and Hernan Arango for his tireless maintenance of the ROMS-IS4DVAR code. The ECCO-GODAE data was provided by the ECCO Consortium for Estimating the Circulation and Climate of the Ocean funded by NOPP. References Baltz, K.A., Sakuma, K.M., Ralston, S., 2006. The physical oceanography off the central California coast during May–June, 2001: a summary of CTD and other hydrographic data from young-of-the-year juvenile rockfish surveys. U.S. Dep. Commer., NOAA Tech. Memo. NMFS, NOAA-TM-NMFS-SWFSC-395, 83 pp. Batchelder, H.P., Barth, J.A., Kosro, P.M., Strub, P.T., Brodeur, R.D., Peterson, W.T., Tynan, C.T., Ohman, M.D., Botsford, L.W., Powell, T.M., Schwing, F.B., Ainley, D.G., Mackas, D.L., Hickey, B.M., Ramp, S.R., 2002. The GLOBEC Northeast Pacific California Current System Program. Oceanography 15 (2), 36–47. Batteen, M.L., 1997. Wind-forced modeling studies of currents, meanders, and eddies in the California Current system. Journal of Geophysical Research 102, 985–1010. Cervantes, B.T.K., Allen, J.S., 2006. Numerical model simulations of continental shelf flows off northern California. Deep-Sea Research II 53 (25–26), 2956–2984. Chapman, D.C., 1985. Numerical treatment of cross-shelf open boundaries in a barotropic coastal ocean model. Journal of Physical Oceanography 15 (8), 1060–1075. Chereskin, T.K., Trunnell, M., 1996. Correlation scales, objective mapping, and absolute geostrophic flow in the California Current. Journal of Geophysical Research-Oceans 104 (C4), 7695–7714. Cummings, J.A., 2005. Operational Multivariate Ocean Data Assimilation. Quarterly Journal of the Royal Meteorological Society 131, 3583–3604. Derber, J.D., Bouttier, F., 1999. A reformulation of the background error covariance in the ECMWF global data assimilation system. Tellus A 51 (2), 195–221, doi:10.1034/j.1600-0870.1999.t01-2-00003.x. Di Lorenzo, E., 2003. Seasonal dynamics of the surface circulation in the Southern California Current System. Deep-Sea Research II 50, 2371–2388. Di Lorenzo, E., Moore, A., Arango, H., Chua, Cornuelle, B.D., Miller, A.J., Powell, B., Bennett, A., 2007. Weak and strong constraint data assimilation in the inverse Regional Ocean Modeling System (ROMS): development and application for a baroclinic coastal upwelling system. Ocean Modeling 16 (3–4), 160–187. Doyle, J.D., Jiang, Q., Chao, Y., Farrara, J., 2009. High-resolution atmospheric modeling over the Monterey Bay during AOSN II. Deep Sea Research II, doi:10.1016/j.dsr2.2008.08.009. Fairall, C.W., Bradley, E.F., Rogers, D.P., Edson, J.B., Young, G.S., 1996. Bulk parameterization of air-sea fluxes for TOGA COARE. Journal of Geophysical Research-Oceans 101, 3747–3767. Fisher, M., 1998. Minimization algorithms for variational data assimilation. In: Proceedings of the ECMWF Seminar on Recent Developments in Numerical Methods for Atmospheric Modelling, Reading, England, 7–11 September 1998, pp. 364–385. Flather, R.A., 1976. A tidal model of the north–west European continental shelf. Memoires de la société Royale des Sciences de Liège 6 (10), 141–164. Gan, J., Allen, J.S., 2002a. A modeling study of shelf circulation off northern California in the region of the Coastal Ocean Dynamics Experiment: Response to relaxation of upwelling winds. Journal of Geophysical Research-Oceans 107 (C9), 3123, doi:10.1029/2000JC000768. Gan, J., Allen, J.S., 2002b. A modeling study of shelf circulation off northern California in the region of the Coastal Ocean Dynamics Experiment 2. Simulations and comparisons with observations. Journal of Geophysical Research-Oceans 107 (C11), 3184, doi:10.1029/2001JC001190. Haidvogel, D.B., Arango, H., Budgell, W.P., Cornuelle, B.D., Curchitser, E., Di Lorenzo, E., Fennel, K., Geyer, W.R., Hermann, A.J., Lanerolle, L., Levin, J., McWilliams, J.C., Miller, A.J., Moore, A.M., Powell, T.M., Shchepetkin, A.F., Sherwood, C.R., Signell, R.P., Warner, J.C., Wilkin, J., 2008. Ocean forecasting in terrain-following coordinates: Formulation and skill assessment of the Regional Ocean Modeling System. Journal of Computational Physics 227, 3595–3624. Heimbach, P., Ponte, R.M., Evangelinos, C., Forget, G., Mazloff, M., Menemenlis, D., Vinogradov, S., Wunsch, C., 2006. Combining altimetric and all other data with a general circulation model. In: Proceedings of the 15 Years of Progress in Radar Altimetry Symposium, Venice, 13–18 March 2006. ESA Special Publication SP-614, ISBN 92-9092-925-1. ESA Publications Division, ESTEC, 2200 AG Noordwijk, The Netherlands. Hickey, B.M., 1979. The California Current System—hypotheses and facts. In: Progress in Oceanography, vol. 8. Pergamon, pp. 191–279. Hickey, B.M., 1998. Coastal oceanography of western North America from the tip of Baja California to Vancouver Island: Coastal segment. In: Robinson, A.R., Brink, K.H. (Eds.), In: The Sea, vol. 11. John Wiley and Sons, pp. 345–391. Hodur, R.M., Pullen, J., Cummings, J., Hong, X., Doyle, J.D., Martin, P., Rennick, M.A., 2002. The coupled ocean/atmosphere mesoscale prediction system (COAMPS). Oceanography 15, 88–98. Kurapov, A.L., Allen, J.S., Egbert, G.D., Miller, R.N., Kosro, P.M., Levine, M., Boyd, T., 2005a. Distant effect of assimilation of moored currents into a model of coastal wind-driven circulation off Oregon. Journal of Geophysical Research-Oceans 110, C02022, doi:10.1029/2003JC002195. Kurapov, A.L., Allen, J.S., Egbert, G.D., Miller, R.N., Kosro, P.M., Levine, M., Boyd, T., Barth, J.A., Moum, J., 2005b. Assimilation of moored velocity data in a model of coastal wind-driven circulation off Oregon: multivariate capabilities. Journal of Geophysical Research-Oceans 110 (c10), C10S08.

92

G. Broquet et al. / Dynamics of Atmospheres and Oceans 48 (2009) 69–92

Large, W.G., McWilliams, J.C., Doney, S.C., 1994. Oceanic vertical mixing: a review and a model with nonlocal boundary layer parameterisation. Reviews of Geophysics 32, 363–403. Le Dimet, F.X., Talagrand, O., 1986. Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects. Tellus 38A, 97–110. Lermusiaux, P.F.J., 2007. Adaptive modeling, adaptive data assimilation and adaptive sampling. Physica D 230, 172–196 (Journal of Inverse Problems, Special issue on “Mathematical Issues and Challenges in Data Assimilation for Geophysical Systems: Interdisciplinary Perspectives”, C.K.R.T. Jones and K. Ide, Eds.). Li, Z., Chao, Y., McWilliams, J.C., Ide, K., 2008a. A three-dimensional variational data assimilation scheme for the Regional Ocean Modeling System. Journal of Atmospheric and Oceanic Technology 25, 2074–2090, doi:10.1175/2008JTECHO594.1. Li, Z., Chao, Y., McWilliams, J.C., Ide, K., 2008b. A three-dimensional variational data assimilation scheme for the Regional Ocean Modeling System: implementation and basic experiments. Journal of Geophysical Research-Oceans 113, C05002, doi:10.1029/2006JC004042. McCreary, J.P.J., Kundu, P.K., Chao, S.-Y., 1987. On the dynamics of the California Current System. Journal of Marine Research 45, 1–32. Marchesiello, P., McWilliams, J.C., Shchepetkin, A., 2001. Open boundary conditions for long-term integration of regional ocean models. Ocean Modelling 3, 1–20. Marchesiello, P., McWilliams, J.C., Shchepetkin, A., 2003. Equilibrium structure and dynamics of the California Current System. Journal of Physical Oceanography 33 (4), 753–783. Moore, A.M., Arango, H.G., Di Lorenzo, E., Cornuelle, B.D., Miller, A.J., Neilson, D.J., 2004. A comprehensive ocean prediction and analysis system based on the tangent linear and adjoint of a regional ocean model. Ocean Modelling 7, 227–258. Moore, A.M., Arango, H.G., Di Lorenzo, E., Cornuelle, B.D., Miller, A.J., 2009. An adjoint sensitivity analysis of the southern California Current circulation and ecosystem. Journal of Physical Oceanography 39, 702–720. Oke, P.R., Allen, J.S., Miller, R.N., Egbert, G.D., Kosro, P.M., 2002. Assimilation of surface velocity data into a primitive equation coastal ocean model. Journal of Geophysical Research-Oceans, doi:10.1029/2000JC00511. Powell, B.S., Arango, H.G., Moore, A.M., Di Lorenzo, E., Milliff, R.F., Foley, D., 2008. 4DVAR Data Assimilation in the Intra-Americas Sea with the Regional Ocean Modeling System (ROMS). Ocean Modelling 25, 173–188. Ricci, S., Weaver, A.T., Vialard, J., Rogel, P., 2005. Incorporating state-dependent temperature-salinity constraints in the background-error covariance of variational ocean data assimilation. Monthly Weather Review 133, 317–338. Shchepetkin, A.F., McWilliams, J.C., 2003. A method for computing horizontal pressure-gradient force in an oceanic model with a nonaligned vertical coordinate. Journal of Geophysical Research-Oceans 108 (C3), 3090, doi:10.1029/2001JC001047. Shchepetkin, A.F., McWilliams, J.C., 2005. The Regional Oceanic Modeling System: a split-explicit, free-surface, topographyfollowing-coordinate ocean model. Ocean Modeling 9, 347–404. Shulman, I., Wu, C.R., Lewis, J.K., Paduan, J.D., Rosenfeld, L.K., Kindle, J.C., Ramp, S.R., Collins, C.A., 2002. High resolution modeling and data assimilation in the Monterey Bay area. Continental Shelf Research 22 (8), 1129–1151. Shulman, I., Kindle, J., Martin, P., deRada, S., Doyle, J., Penta, B., Anderson, S., Chavez, F., Paduan, J., Ramp, S., 2007. Modeling of upwelling/relaxation events with the Navy Coastal Ocean Model. Journal of Geophysical Research-Oceans 112, C06, 023, doi:10.1029/2006JC003946. Taylor, K.E., 2001. Summarizing multiple aspects of model performance in a single diagram. Journal of Geophysical Research 106 (D7), 7183–7192. Tshimanga, J., Gratton, S., Weaver, A.T., Sartenaer, A., 2008. Limited-memory preconditioners with application to incremental variational data assimilation. Quarterly Journal of the Royal Meteorological Society 134, 751–769. Veneziani, M., Edwards, C.A., Doyle, J.D., Foley, D., 2009a. A Central California coastal ocean modeling study: 1. Forward model and the influence of realistic versus climatological forcing. Journal of Geophysical Research-Oceans 114, C04015, doi:10.1029/2008JC004774. Veneziani, M., Edwards, C.A., Moore, A.M., in press-b. A Central California modeling study. Part II: Adjoint sensitivities to local and remote driving mechanisms. Journal of Geophysical Research-Oceans, doi:10.1029/2008JC004775. Warner, J.C., Sherwood, C.R., Arango, H.G., Signell, R.P., 2005. Performance of four turbulence closure models implemented using a generic length scale method. Ocean Modeling 8, 81–113. Weaver, A.T., Courtier, P., 2001. Correlation modelling on the sphere using a generalized diffusion equation. Quarterly Journal of the Royal Meteorological Society 127, 1815–1846. Weaver, A.T., Vialard, J., Anderson, D.L.T., 2003. Three- and four-dimensional variational assimilation with a general circulation model of the tropical pacific ocean. Part I: Formulation, internal diagnostics, and consistency checks. Monthly Weather Review 131, 1360–1378. Willis, J.K., Roemmich, D., Cornuelle, B., 2004. Interannual variability in upper ocean heat content, temperature, and thermosteric expansion on global scales. Journal of Geophysical Research-Oceans 109, C12036, doi:10.1029/2003JC002260.