Weak constraint four-dimensional variational data assimilation in a model of the California Current System

A new approach is explored for computing estimates of the error covariance associated with the intrinsic errors of a numerical forecast model in regions characterized by upwelling and downwelling. The approach used is based on a combination of strong constraint data assimilation, twin model experiments, linear inverse modeling, and Bayesian hierarchical modeling. The resulting model error covariance estimates Q are applied to a model of the California Current System using weak constraint four-dimensional variational (4D-Var) data assimilation to compute estimates of the ocean circulation. The results of this study show that the estimates of Q derived following our approach lead to demonstrable improvements in the model circulation estimates and isolate regions where model errors are likely to be important and that have been independently identified in the same model in previously published work.


Introduction
Data assimilation has been a mainstay of numerical weather prediction since the 1980s and is a critical component of any forecast system.Data assimilation methods are also used routinely for computing estimates of the ocean circulation based on sparse and incomplete observations, and have been adopted in other fields such as sea-ice modeling and biogeochemical modeling of the atmosphere and ocean.A thorough review of data assimilation methods and applications in the geosciences is beyond the scope of this paper, but excellent reviews can be found in Daley (1991), Bennett (2002), Kalnay (2002), and Wunsch (2006).
Formally, data assimilation methods can be derived from Bayes' theorem (Wikle and Berliner, 2007) and can be viewed as the combination of a prior estimate (usually from a forecast model) and observations, appropriately weighted according to a priori assumptions about the relative uncertainties in both.The Kalman filter (or smoother) provides the most complete solution to linear data assimilation problems.
However, the dimension of most atmosphere and ocean models prohibits application of the Kalman filter in its full form, and some approximations or alternative approaches must be considered.Data assimilation methods that are now in common use at operational forecast centers for computing atmospheric and oceanic analyses generally fall into three categories: (i) ensemble Kalman filters, (ii) variational methods, and (iii) hybrid combinations of (i) and (ii).Ensemble Kalman filters employ Monte Carlo methods to estimate (and propagate in time) the prior error and posterior error covariances employed in the full Kalman filter formulation.Variational methods, on the other hand, employ the methods of optimal control theory to directly identify the state of the system that corresponds to the maximum of the posterior probability distribution.Both approaches have important strengths (and weaknesses) and hybrid methods are an attempt to capitalize on these.In all cases, however, a linear analysis is assumed subject to errors that belong to a Gaussian distribution.Our focus here is on variational data assimilation methods since they are employed in the present work.
Regardless of the variety of data assimilation method employed, most operational centers ignore the direct influence of model errors when computing an analysis.This assumption is made not because the forecast models are believed to be perfect, but because the nature and size of model errors is very poorly understood, which hampers any formal treatment of them during the analysis procedure.In this study, we demonstrate a dynamically based method for estimating a reduced-rank form of the model error covariance of an ocean model that can be used during data assimilation in regions where the ocean circulation and wind forcing are tightly coupled.
The paper is structured as follows.In Sect.2, we present a brief overview of the variational data assimilation method used in our experiments.The ocean model configuration is described in Sect.3, while Sect. 4 describes in some detail the methods used to estimate the model error covariance matrix.The resulting estimates of the model error covariance are tested via a series of weak constraint variational data assimilation experiments which are described in Sect. 5. A summary and discussion of the findings from the experiments follows in Sect.6.

4D-Var data assimilation
In this study, we employ a four-dimensional variational (4D-Var) data assimilation method for estimating the timeevolving state of the ocean circulation.Following the standard notation used in the meteorological and oceanographic literature (Ide et al., 1997;Daget et al., 2009), the ocean state vector will be denoted as x and comprises the model grid point values of temperature, salinity, two components of velocity, and sea surface height.The state vector x(t i ) is advanced forward in time by the forecast model, which is denoted as M so that where x(t i ) represents the initial state, and f (t i , t f ) and b(t i , t f ) represent the ocean surface forcing and lateral boundary conditions, respectively, for the time interval t = [t i , t f ].All observations during the same interval will be denoted as the elements of the vector y.If x b (t i ), f b (t i , t f ), and b b (t i , t f ) denote the background estimates of the circulation, surface forcing, and lateral boundary conditions, the goal of 4D-Var is to identify the analysis estimates x a (t i ), f a (t i , t f ), and b a (t i , t f ) that minimize the cost function where z = (x T (t i ), f T (t i , t f ), b T (t i , t f )) T is the vector of control variables, and H denotes the observation operator that maps z into the observation space.In general, H will also include the model M. Uncertainties in the background control vector z b (t i ) and the observations y are described by the background error and observation error covariance matrices B and R, respectively.For Gaussian errors, the conditional probability of z given z b and y is proportional to e −J NL , in which case the z that minimizes J NL also maximizes the conditional probability (Wikle and Berliner, 2007).The calculus of variations can be used to identify the z that minimizes J NL (Le Dimet and Talagrand, 1986).However, because the nonlinear forecast model M is an integral component in H in 4D-Var, the cost function J NL may have a complicated topology which can make identification of the minimum difficult.It is therefore common practice to linearize the state estimation problem about the background under the assumption that increments δz to z b are small, so that z a = z b + δz.Following this so-called incremental approach (Courtier et al., 1994), the variational data assimilation problem is transformed to one in which an iterative sequence of linear minimizations are performed to identify the increments δz k , where the cost function is given by where k is iteration index, d k−1 = y − H (z k−1 ) is referred to as the innovation vector, and G k−1 is the generalized observation operator and represents a tangent linearization of H about z k−1 (t).The background error covariance matrix B k in general depends on the outer-loop iteration via the presence of a state-dependent balance operator (see Sect. 3).In 4D-Var, G k−1 represents the tangent linear model sampled at the observation locations and G T k−1 the adjoint model forced at the observation locations.The iterates z k are referred to as outer loops, and each iterate is identified by minimizing J k via a sequence of inner loops.After each outer loop, a new iterate is computed according to z k = z k−1 + δz k .The final iterate in the sequence identifies the analysis circulation estimate z a = z b + n k=1 δz k .The iterative procedure involving inner and outer loops is equivalent to a truncated Gauss-Newton method (Lawless et al., 2005) and yields an estimate of the control vector z that minimizes J NL .
As currently formulated, the minimization of J NL in Eq. ( 2) via the sequence of inner and outer loops of the form Eq. ( 3) neglects errors in the model M. As a result, the resulting analysis circulation estimate x a is constrained to be an exact solution of the model equations.This is commonly referred to as strong constraint data assimilation (Sasaki, 1970) and is the practice frequently adopted at most operational centers.However, the control vector z can be augmented to include the influence of model errors also, so that δz = (δx T (t i ), δf T (t i , t f ), δb T (t i , t f ), η T (t i , t f )) T , where η(t i , t f ) represents corrections to the model grid point values over the interval t = [t i , t f ] to account for the presence of model error.In this case, the background error covariance matrix B now contains a block Q along the leading diagonal that describes the covariance of the model errors.Specifically, B = diag(B x , B f , B b , Q), where B x , B f , and B b are the initial condition, surface forcing, and lateral boundary condition background error covariance matrices.The block diagonal structure of B assumes that the errors in each control variable type are uncorrelated (i.e., the initial condition errors are assumed to be uncorrelated with the surface forcing errors), although each block will typically be a non-diagonal, fullrank covariance matrix.This approach differs fundamentally from strong constraint 4D-Var in that the resulting analysis circulation estimate x a is now only an approximate solution of the model equations M. The dynamics of the model are therefore imposed on the estimate as a weak constraint (Bennett, 2002).
Both strong and weak constraint data assimilation present considerable challenges in the large dimensional systems encountered in atmosphere and ocean modeling.The curse of dimension, however, is exacerbated by the weak constraint formulation because of the need to compute the elements of η at every model grid point and possibly every time step.This problem can be mitigated by using the dual approach to 4D-Var.The iterate z k that minimizes J k in Eq. ( 3) is given by where K k−1 is the Kalman gain matrix for the kth outer loop.The Kalman gain matrix can be written in two equivalent forms: or Equation ( 5) is referred to as the primal form and involves computing an estimate of the inverse of the Hessian matrix which has the dimension of the control space.Conversely, Eq. ( 6) is referred to as the dual form and involves computing an estimate of the inverse of the stabilized representer matrix (G k B k G T k + R) which has a dimension equal to the number of observations, regardless of the size of the control vector.Therefore, the dual form presents a more tractable approach to the weak constraint 4D-Var problem and will be used here.

Model configuration
The model used here is the Regional Ocean Modeling System (ROMS), which solves the primitive equations subject to the hydrostatic approximation.ROMS is designed primarily for regional applications (Haidvogel et al., 2008) and employs curvilinear coordinates in the horizontal direction and terrain-following coordinates in the vertical direction, making it very flexible for coastal applications (Shchepetkin and McWilliams, 2005).A state-of-the-art 4D-Var system is also available for ROMS, both in the primal and dual formulations (Moore et al., 2011b).Using the same background estimates, observations, and error covariances, both the primal and dual forms of ROMS 4D-Var yield the same circulation estimates (Gürol et al., 2013).
The configuration of ROMS employed here is shown in Fig. 1 and spans the California Current System (CCS) which is the dominant circulation environment off the west coast of North America.The CCS forms the eastern branch of the North Pacific subtropical gyre and is subject to pronounced seasonal variations (Hickey, 1998).During the spring and summer, the winds along the entire coast from the Canadian border to Baja California are equatorward.This drives offshore Ekman transport in the surface layers of the ocean creating cold, nutrient-rich conditions immediately adjacent to the coast.These conditions are favorable for blooms of phytoplankton and fuel high levels of primary production in the ocean, which in turn support a variety of rich marine ecosystems (Checkley and Barth, 2009).A schematic of some important circulation features in the region is shown in Fig. 1.In fact, the CCS is one of four so-called eastern boundary large marine ecosystems that also include the Humboldt-Peru Current in the South Pacific, the Canary Current in the North Atlantic, and the Benguela Current in the South Atlantic.As such, ocean state estimation and forecasting of the CCS is of considerable socioeconomic importance.
The model was configured with 1/10 • horizontal resolution and 42 terrain-following σ levels in the vertical that vary in thickness between 0.3 and 8 m over the continental shelf and between 7 and 100 m in the deep ocean.
ROMS strong constraint 4D-Var has been successfully applied to the CCS in a series of studies (Broquet et al., 2009a(Broquet et al., , b, 2011;;Moore et al., 2011aMoore et al., , 2013) ) culminating in two long sequences of historical circulation analyses spanning the period 1980-2012 (Neveu et al., 2016, hereafter N16).Experience with weak constraint 4D-Var is more limited (Moore et al., 2011a) because of the difficulty in quantifying model error.
The present study capitalizes on the WCRA14 historical analyses of N16.Therefore, only a brief description of the 4D-Var configuration will be given here, and interested readers should consult N16 and references therein for more detailed information.The specific N16 analysis period considered here is 1999-2012.During this time, the background surface forcing f b was derived from near-surface atmospheric fields from the Coupled Ocean-Atmosphere Mesoscale Prediction System (COAMPS) described by (Doyle et al., 2009).As shown in Fig. 1  The observations y take the form of satellite-derived sea surface temperatures (SSTs) from the AVHRR, MODIS and AMSR-E instruments, a satellite-derived gridded sea surface height product from Aviso (specifically, the DUACS version DT-2010), and quality-controlled in situ profiles of temperature and salinity (confined primarily to the upper 1000 m) from the Met Office EN3 database (v2a) described by Ingleby and Huddleston (2007).Observation errors were assumed to be uncorrelated, in which case R is a diagonal matrix.While this is a reasonable assumption for in situ measurements, it is unlikely to be correct for satellite observations.However, incorporating the effects of correlated observation errors in current data assimilation systems is very challenging and currently a topic of active research.
Following Weaver et al. (2005), the circulation increments δx can be decomposed into dynamically balanced and unbalanced components, and the background error covariance matrix for the increments factorized according to The unbalanced increments are assumed to have no cross-correlations and are described by the univariate correlation matrix C, and the standard deviations are represented by the elements of the diagonal matrix .Cross covariances are introduced by the balance operator K b .The univariate correlation matrix C is modeled as the solution of a pseudo-heat diffusion operator following Weaver and Courtier (2001).The background error covariances B f and B b are factorized in a similar way, except that there is no balance operator.A discussion of the choice of correlation length scales and standard deviations used to model each of the background error covariances can be found in N16.The balance operator was not used in N16.
Strong constraint 4D-Var was applied sequentially by N16 over the 1999-2012 period using overlapping 8-day assimilation windows.All available observations during each window were assimilated into the model.Each assimilation window overlaps with the previous window by 4 days, so the analysis circulation estimate x a at the midpoint of the previous window was used as the background circulation x b for the current assimilation window.In all cases, 1 outer loop and 15 inner loops were used to approximate the minimum of J NL using Eq.(3) and the dual formulation of 4D-Var (i.e., circulation estimates based on Eq. 4 using Eq. 6).
Using a non-data-assimilative configuration of the same model, Veneziani et al. (2009) showed that the difference between modeled and observed SSTs along the central and northern California coast varies seasonally.During the winter and spring, the model SST between Cape Mendocino and Point Conception (Fig. 1) is typically colder than observed by ∼ 0.5 • C, although in spring the model is warmer by ∼ 0.25 • C right at the coast south of Point Arena.Meanwhile, during the fall the entire central California SST is warmer than observed.North of Cape Mendocino, the model is typically colder than observed year-round by ∼ 0.5-1 • C. The standard deviation of the model minus observation differences in SST were found to be smallest during the spring and largest during fall, when the eddy kinetic energy in the region also reaches a maximum.During strong constraint 4D-Var experiments in which the forcing was included in the control vector, Broquet et al. (2011) found that a reduction in the model minus observed difference in SST was accompanied by a change in the strength of the alongshore surface wind stress derived from COAMPS.Changes in the strength of alongshore winds during 4D-Var led to changes in the rate of upwelling in the model and surface temperatures in closer agreement with observations, although changes in the 4D-Var adjusted surface heat flux were also found to be important in some areas.The COAMPS wind stress agrees very well with satellite scatterometer observations in the region (Doyle et al., 2009), while the 4D-Var corrected wind stress is typically weaker than observed compared to the same scatterometer product (Broquet et al., 2011).The associated changes in surface heat flux are of course harder to quantify due to the lack of good observations.However, based on these previous studies, it appears that during data assimilation the influence of errors in the model on the upwelling circulation are compensated for by changes in the surface forcing.It is these findings that motivate the present study and the methodology that is proposed in Sect. 4. Similar results have been reported by Stammer et al. (2002), who found that the addition of vorticity via corrections to wind stress in the vicinity of boundary currents in a coarse-resolution ocean model is needed to yield proper separation of currents from the coast.
Much of the aforementioned primary coastal upwelling circulation along the North American west coast can be understood in terms of linear dynamics (e.g., Gill, 1982).While the rate of upwelling is controlled by the coastal divergence of Ekman transport, the resulting vertical displacement of isopycnals depends also on the stratification (Gill, 1982), and interactions with the bottom boundary over the continental shelf can further complicate the circulation (Jacox and Edwards, 2011).However, if we view coastal upwelling in the CCS as a predominantly linear process, then we argue (based on the previous aforementioned findings) that errors in the ocean model temperatures in coastal upwelling zones that are associated with errors in the model formulation (e.g., numerics and parameterizations) may be compensated for during 4D-Var by adjustments in the surface wind stress.Linear dynamics would seem to imply a one-to-one relationship between the surface wind stress and coastal ocean response.So, conversely, if the ocean model is subjected to the 4D-Var adjusted wind stress, the upwelling response will be affected in such a way as to mitigate the temperature errors associated with model error.This suggests that in a linear dynamical regime, the 4D-Var corrections made to good quality ocean forcing estimates f should provide direct information about the circulation corrections that are required to offset the influence of intrinsic model errors.In practice, because the two-dimensional surface forcing f projects onto a three-dimensional circulation field x, a unique one-to-one relationship will generally not exist between the two.In the "forward" problem, the surface forcing f does uniquely determine the 3-D circulation x (in conjunction with the initial conditions and boundary conditions, of course).However, in the "inverse" problem, given the 3-D circulation x, the surface forcing f cannot, in general, be uniquely determined unless additional constraints are imposed on the problem.This is the strategy employed by data assimilation methods such as 4D-Var.Specifically, 4D-Var identifies corrections to the prior surface forcing f b (t) subject to the constraints imposed by prior information (i.e., the model dynamics M, the initial conditions x b (t i ), boundary conditions b b (t i , t f ), background error B, and observation error R covariance matrices).If we have confidence in the prior information, then we can place equal confidence in the 4D-Var posterior estimates of the surface forcing.Therefore, 4D-Var alleviates the underdetermined nature of the inverse relationship between the surface forcing and the coastal upwelling response, and it is this property that we aim to capitalize on here.

Estimates of the model error covariance
Forecast models are subject to model errors (e.g., discretization errors, errors associated with parameterizations of unresolved processes), and as noted in Sect.2, these errors should be accounted for during data assimilation.However, the nature and magnitude of model errors is poorly known, making it difficult to reliably estimate the model error covariance matrix Q.In some cases, careful consideration of the various sources of errors can be fruitful (e.g., Bennett, 2002, and references therein).The method used here to estimate sources of model error differs from previous approaches and is based on a combination of twin model experiments (Sect.4.1), linear inverse modeling (Sect.4.2), and Bayesian hierarchical modeling (Sect. 4.3).
By way of an introduction to the approach used here, consider a perfect ocean forecast model subject to error-free surface forcing, so that where N represents the model operators, w t (t) is the influence of the surface forcing on the ocean interior (e.g., the projection of the surface forcing onto the barotropic and baroclinic modes), and x t denotes the true ocean circulation.Specifically, w(t) is the result of a linear map of the surface forcing f onto x via the boundary conditions for the vertical diffusion of momentum, heat, and salt.For example, for the zonal momentum, the surface boundary conditions are κ∂u/∂z = (1/ρ o )τ x , where τ x is the zonal surface wind stress and κ is the vertical diffusion coefficient.Therefore, we can represent w(t) as w(t) = Pf (t), where P is the projection matrix that depends on the vertical mixing scheme that is employed.We will assume that P is independent of the state x; however, P could be a function of x depending on the vertical mixing parameterization used.Consider now an imperfect model, but still subject to error-free surface forcing: where (t) represents the error in the time tendency of each element of the state vector.The assumption here is that (t) takes the form of additive noise and is independent of the state vector x.In general, this may not be true of errors in model parameterizations, but it is a reasonable starting point here.However, the more challenging problem of multiplicative noise deserves further attention.From Eqs. ( 7) and ( 8), the evolution of the state-vector differences x = x−x t over time will be given by where the second equality follows from a first-order Taylor expansion and N = ∂N /∂x| x t .When N is autonomous1 , solutions of Eq. ( 9) can be written as which shows that the state-vector differences at any time depend on the time-evolved initial error x(0) and the integrated influence of the model error, (t).
Consider the following thought experiment, and suppose that we perform strong constraint 4D-Var using the imperfect forecast model Eq. ( 8), with error-free surface forcing, and using the initial conditions and surface forcing as the control variables, z.The data assimilation procedure will attempt to compensate for the presence of model errors through erroneous adjustments to both the initial conditions and surface forcing.
Consider now a pair of model integrations using Eq. ( 8) (a "twin model" approach), one subject to error-free forcing, w t , that yields the circulation x 1 , and a second subject to the 4D-Var corrected forcing, w a , that yields the circulation x 2 .To first order, the circulation difference δx = x 2 − x 1 will be given by where M = ∂N /∂x| x 1 and δw = w a − w t .Equation ( 11) clearly has the same mathematical form as Eq. ( 9).Therefore, if the two model integrations for x 1 and x 2 start from the same initial condition (x 1 (0) = x 2 (0)), we hypothesize that the circulation differences δx that develop over time in response to δw = P(f a − f b ) will mimic the characteristics and influence that the model errors (that give rise to δw during strong constraint 4D-Var) have on the circulation error x.That is to say, there will be some model error (t) which will yield the circulation errors x(t), and it is this (t) which we aim to identify.In general, the model error will also project onto the initial conditions (see Appendix A), but the approach we propose here will identify only the signature of model error in the surface forcing.As described in Sect.4.3, a critical element of the approach adopted here is the partitioning of the increments δz between legitimate corrections and those associated with model error.This requires independent estimates of the probability distribution of each component of the control vector.Such estimates are available for f using Bayesian methods (see Sect. 4.3), but not for x(t i ) due to the large dimension and multivariate nature of the ocean state.In fact, the goal of data assimilation is to estimate properties of the state-vector distribution.Therefore, we confine our attention here to what we can learn about model errors from the surface forcing increments.
To estimate the model error covariance E{ T } = Q, we use the technique of linear inverse modeling as described in Sect.4.2.While the surface forcing will never be known precisely, uncertainties in the forcing fields can be accounted for if their probability distributions can be estimated.To this end, we will use the methods of Bayesian hierarchical modeling as described in Sect.4.3.

Twin model approach
The "twin model" approach outlined above was applied to ROMS CCS circulation estimates from the WCRA14 historical analyses described by N16 using strong constraint 4D-Var.As noted in Sect.3, these analyses cover the period 1999-2012, and each 4D-Var assimilation window spans an 8-day interval.In the following, we confine our attention to the shorter period [2003][2004].Two sequences of nonlinear ROMS integrations were performed spanning this period.During the first sequence, the model was initialized with the 4D-Var background circulation for each 8-day assimilation window and integrated forward in time for 8 days subject to the COAMPS background surface forcing.The circulations from this sequence will be denoted by x i , where i denotes the 4D-Var cycle number.During the second sequence, the model was initialized with the same background circulation as in x i , but in this case the strong constraint 4D-Var analysis estimates of the wind forcing were used to integrate the model forward in time over the same 8-day window.These winds represent the corrections made by strong constraint 4D-Var to the COAMPS background and, as argued earlier, contain a signature of the model error (see also Appendix A).The background fluxes of heat and freshwater were constrained to be the same in both sequences.The circulations from the second sequence will be denoted by xi .In this demonstration study, the focus is on the surface wind forcing alone because only a Bayesian hierarchical model for the winds is available at the present time (see Sect. 4.3).The circulation sequences x i and xi are synonymous with x 1 and x 2 described earlier in connection with Eq. ( 11).During each 4D-Var cycle, an estimate of the ocean circulation is computed every 8 days as described by N16.The corrections that must be made to account for the presence of model error during each 4D-Var cycle should therefore reflect the growth and evolution of circulation changes due to model errors during an 8-day period, rather than the climatological variance associated with such errors.For example, in the case of perfect initial conditions, surface forcing, and open boundary conditions, the errors in the circulation on day 8 will be those associated with model error only during that period.This is the rationale for focusing on the circulation differences on day 8 in the twin experiments.

Linear inverse modeling and estimation of Q
Assume for the moment that the surface fluxes derived from the COAMPS near-surface fields are error-free.They are not, of course, but we will relax this assumption in Sect.4.3.By the aforementioned logic, the sequence of model circulation differences δx i = xi − x i from Sect.4.1 can be viewed as an expression of the influence of model error, (t), on the CCS circulation estimates.In this sense, we can view the δx i as an ensemble of the influence of different realizations of the model error (t) on the circulation.Here, we focus our attention on the δx i at the end of each 8-day assimilation window since, as noted above, we are interested in estimating the influence of model errors over the length of a typical assimilation window.By invoking the ergodic hypothesis, we can view the ensemble δx i as equivalent to a time series δx(t i ) and make a further assumption that this time series can be modeled as a first-order Markov process.Specifically, we will assume that the δx vectors are described by where A is a linear operator that advances δx forward one time step dt = t i+1 − t i , and ξ (t i ) is stochastic forcing that is white in time.Equation ( 12) is a discrete analog of Eq. ( 9) and reveals that, apart from a scaling factor, the stochastic forcing ξ (t i ) is synonymous with the model error (t).While it is unlikely that the true model error (t) is an additive white noise time-continuous process, this may be a reasonable assumption for the discrete case considered here depending on the choice of dt, which in the current case is 4 days (i.e., the time between each overlapping 8-day assimilation window).
We will return to this point shortly.Following von Storch et al. (1995), right multiplying both sides of Eq. ( 12) by δx T (t i ) and taking the expected value, yields ) where C 0 and C 1 are the lag-0 and lag-1 covariance matrices of δx(t i ), respectively.Once A has been estimated from Eq. ( 13), the stochastic forcing ξ can be diagnosed from Eq. ( 12) according to ξ (t i ) = (δx(t i+1 ) − Aδx(t i )).Formally, the model errors (t) in Eq. ( 9) can be represented as a Wiener process (i.e., Brownian motion).For the discrete form of Eq. ( 9) that is equivalent to Eq. ( 12), (t (Kloeden and Platen, 1992), where dt = t i+1 − t i is assumed constant (i.e., the mean squared distance between two realizations of the Wiener process increases linearly with time).Therefore, scales as ξ /dt 1 2 , where it is understood that the square root is applied to the numerical value of dt and not the time units.Following Penland and Matrosova (1994), the model error covariance matrix can then be estimated according to Q = E{ξ ξ T }/dt.As noted in Sect. 1, the dimension of the state vector of atmospheric and oceanic forecast models is typically very large and, in the case of ROMS CCS, is ∼ 10 6 − 10 7 .This precludes the direct computation of A based on Eq. ( 13).However, as shown by von Storch et al. (1995) and Penland and Sardeshmukh (1995) the dimension of the problem can be greatly reduced by expanding δx in terms of the empirical orthogonal functions (EOFs) of C 0 .Following Crawford et al. ( 2016), Eq. ( 12) was transformed into an equivalent first-order Markov equation for the variable u = W 1 2 δx, where W is the matrix such that δx T Wδx defines the physical energy per unit volume (hereafter referred to as the energy density) of δx.In this way, all of the state variables are appropriately scaled and the trace of the lag-0 covariance matrix C 0 = E{u(t i )u T (t i )} yields the total energy density.This choice of W was predicated on the fact that energy is a fundamental physical quantity that naturally links all elements of the state vector.
The leading EOFs and eigenvalues of C 0 were computed using a Lanczos algorithm (Golub and van Loan, 1989) as described by Crawford et al. (2016) which negates the need to explicitly compute the matrix C 0 .The Lanczos algorithm is iterative, and Fig. 2a shows the leading eigenvalue estimates of C 0 that result from 120 iterations.Figure 2b shows a formal estimate of the error | C 0 ψ i − λ i ψ i | in each estimated eigenpair (λ i , ψ i ).The leading 40 or so eigenpairs have very small expected errors (< 10 −10 ), and even for EOF 50 we consider the errors to be acceptably small.Figure 2c shows the cumulative percentage variance explained by the EOFs.In the following analyses, the leading N = 50 EOFs were used to compute estimates of Q and account for ∼ 61 % of the energy density of δx i .Spatial smoothing was used to suppress any spurious covariances that may result from a finitesized ensemble.If the δx values are spatially smoothed using five applications of a second-order Shapiro filter prior to computation of C 0 , the variance explained by the leading 50 EOFs increases to ∼ 75 % (Fig. 2c).
If denotes the matrix of the leading N EOFs, then δx(t i ) W − 1 2 p(t i ), where p(t i ) is the vector of N principal component time series for each EOF.Equation ( 12) can then be recast in terms of p(t i ) so that   where Â and ξ are the system matrix and stochastic forcing in the reduced dimensional space.Analogous to Eq. ( 13), Â can be estimated according to where 1 and 0 are the lag-1 and lag-0 covariance matrices for the principal component time series p(t).
Figure 3 shows the autocorrelation for each component of ξ for N = 50 using unsmoothed δx i .In each case, the stochastic forcing decorrelates very quickly in time and within one time step, dt, as expected for a Wiener process.We conclude therefore that Eq. ( 13) is an appropriate model for δx.The model error covariance matrix is approximated as , where Q = E{ ξ ξ T }/dt.Because of the linear relationship between coastal upwelling/downwelling and alongshore wind stress, the ap- The eigenspectrum estimate of Q and the cumulative variance explained by each eigenvector is shown in Fig. 6 for smoothed and unsmoothed δx.The surface structure of the leading EOF of Q is shown in Fig. 7, explains ∼ 10 % of the variance, and shares many features in common with the total variance of Fig. 4.

Bayesian hierarchical modeling of surface winds
The discussion in preceding sections has progressed based on the assumption that the COAMPS surface forcing is errorfree and perfectly known.This is of course not the case, even given the high degree of correspondence between the COAMPS wind fields and independent observations (Doyle et al., 2009).Therefore, before performing the twin model experiments of Sect.4.1, it is desirable to obtain an estimate of the uncertainty in the COAMPS background surface forc-ing.Based on these estimates of the uncertainty in the background forcing, the analysis surface forcing estimates from 4D-Var can then be assessed to determine if they are consistent with the surface wind distribution of COAMPS.Specifically, consider the following cases: -Case 1: if the analysis surface forcing estimates from 4D-Var are statistically indistinguishable from a sample drawn from the surface wind distribution of COAMPS, then we accept the analysis estimate as a legitimate representation of the COAMPS state.
-Case 2: if the analysis surface forcing estimates from 4D-Var fall outside the COAMPS probability distribution, then we would reject this as a likely state of COAMPS and attribute the erroneous analysis forcing to the influence of model error.
These comparisons can be performed at each model grid point.In terms of the twin model run for x 2 , at grid points that satisfy Case 1, we accept the background COAMPS forcing as an good estimate of the true forcing, while at grid points that satisfy Case 2, we replace the COAMPS forcing with the 4D-Var analysis estimate of surface forcing.In this way, x 2 (t) will be subject to the background COAMPS surface just as in x 1 except at locations where the COAMPS surface wind distribution indicates that the 4D-Var posterior surface forcing estimates are unduly influenced by model error.In theory, the same approach could be applied to the signature of model error in the initial conditions as well, and as demonstrated in Appendix A, this would be an effective method for recovering model error information over the entire domain.However, as previously noted, estimates of the distribution of the ocean circulation are not available, so in practice there is no way to distinguish between legitimate initial condition corrections to x(t i ) during strong constraint 4D-Var and those associated with the signature of model error.Therefore, as noted earlier, in this demonstration study we will restrict our attention to the surface wind stress components from COAMPS.The surface wind distribution for the COAMPS 4D-Var background wind forcing was estimated using a Bayesian hierarchical model (BHM).The BHM wind model for the CCS is modified slightly from the implementation for the Mediterranean Sea as described in Milliff et al. (2011) (see their Appendix for a complete model specification).The important difference in the BHM design for the CCS involves a change from the multi-resolution nested wavelet bases used to enforce a k −2 slope for the surface wind component of kinetic energy vs. spatial wave number, at the highest resolved wave numbers.In the CCS implementation, the wavelet bases are replaced with Fourier modes.Hyperprior amplitude specifications for a sequence of highwave-number two-dimensional Fourier modes enforce the desired kinetic energy spectral behavior.
The surface wind BHM for the CCS ingests higher-spatialresolution data stage inputs as well.Surface pressure and surface wind components (u, v) are taken from COAMPS anal-yses at 0.1 • resolution, at four canonical times each day, i.e., 03:00, 09:00, 15:00, and 21:00 UTC.Surface vector wind observations are taken from in-swath retrievals in the Level 2B data set for the QuikSCAT mission as provided by the JPL PODAAC (http://podaac.jpl.nasa.gov/dataset/QSCAT_LEVEL_2B_OWV_COMP_12).The QuikSCAT winds are produced at 12.5 km resolution within the observation swath, and QuikSCAT swaths cross the CCS domain twice each day near 03:00 and 15:00 UTC.We assign the QuikSCAT data to these times in the BHM.
Given the COAMPS and QuikSCAT input data stages, the CCS surface wind BHM produces daily estimates of the posterior mean surface wind vector as well as 10 realizations from the posterior distribution at each grid location on a 25 km regular grid from 30 to 42 • N and from 135 • W to the coast.Data stage inputs are weighted such that inputs obtained closest to the output time of 12:00 UTC are twice as influential as earlier and later inputs in the iterations leading to the posterior distribution.
Figure 8 shows the QuikSCAT and COAMPS surface vector wind inputs as assigned to 03:00, 09:00, 15:00, and 21:00 UTC in the CCS domain.The observations from 09:00 and 15:00 UTC are twice as important in the data stage distribution as the observations from 03:00 or 21:00 UTC because they are closer in time to the nominal forecast model output time of 12:00 UTC for 19 January 2003.
The abundance of domain-filling and relatively precise observations from QuikSCAT have a marked effect on surface wind uncertainty estimates.While snapshots of the surface vector wind posterior mean and realizations represent meteorological variability at atmospheric mesoscale and synoptic scales resolved on the 25 km grid, the estimates of the variance show data stage impacts as well, as shown in Fig. 9a  and b.Note the swath patterns in the surface vector wind variance estimates in Fig. 9a and b.
The parameters of the pressure-gradient wind balance imposed in the process model part of the BHM (see Milliff et al., 2011) are estimated as part of the posterior distribution.Traces over the course of the Gibbs iterations (not shown) demonstrate that the BHM has converged to the target posterior distribution, and samples from this portion of the chains are valid samples from the posterior distribution.For example, the 10 realizations shown at each grid location in Fig. 9c are selected after the BHM chains have converged.Importantly, to ensure that these are independent samples, each realization is separated by 10 000 iterations starting after the BHM has converged to the posterior distribution.
In the twin model experiments of Sect.4.1, surface wind stress derived from COAMPS was assumed to be error-free.This assumption will be relaxed here using the CCS BHM for surface winds.Our working hypothesis has been that during strong constraint 4D-Var, corrections to the surface forcing are a manifestation not only of uncertainties in the model forcing but also of errors in the model in regions of coastal upwelling and downwelling.Since the BHM pro-   vides an estimate of the posterior space-time wind distribution, we can use this information to make an informed decision about where the 4D-Var corrected winds are unduly influenced by model error.Specifically, in cases where the 4D-Var corrected winds fall within the BHM posterior distribution, the 4D-Var winds represent a physically realizable state of the atmosphere, and the influence of model error is assumed to be minimal.Conversely, if the 4D-Var corrected winds fall outside the BHM distribution, this is taken to represent a situation where model error exerts a considerable influence on the ocean circulation and where 4D-Var is pushing the winds into a state that is not realizable but necessary in order to better fit ROMS to the observations and circulation background estimate.Following this approach, a second set of twin model experiments was performed spanning the period 2003-2004.The first sequence was identical to x i described in Sect.4.1 in which the model was initialized with the 4D-Var background circulation for each 8-day window of N16 and integrated forward in time for 8 days subject to the COAMPS background surface forcing.During the new second sequence, the model was initialized with the same 4D-Var background circulation for each 8-day window, but in this case forced with surface wind fields that were a combination of the COAMPS background and the 4D-Var adjusted winds, hereafter referred to as "mixed winds".At model grid points and times where the 4D-Var adjusted increments fall within the BHM posterior distribution (within ±2 standard deviations), COAMPS winds were used; otherwise, the 4D-Var adjusted winds were used.For example, Fig. 9c shows a comparison of the 4D-Var adjusted wind stress (black vectors) with 10 realizations of surface stress (vector cluster in red) derived from the surface wind BHM within the subdomain 128-124 • W, 40-44 • N.This is a region where the 4D-Var corrections to the COAMPS background wind stress are generally large and where model error is known to be important (see Fig. 10).The sequence of circulation estimates derived from the mixed winds will be referred to as x i .The surface heat and freshwater fluxes were constrained to be the same during both sequences of runs.
Figure 10 shows the root mean square (rms) difference between the COAMPS background wind stress and the mixed wind stress and reveals that the largest differences typically occur close to the coast near to and equatorward of Cape Blanco.As described in Sect.4.2, the circulation differences δ x i = x i −x i were represented by a first-order Markov model (Eq.12), and a new estimate of the model error covariance matrix Q was computed.The eigenspectra of A and Q are similar to those described in Sect.4.2 (not shown).
It should be noted that the approach used here to compute the mixed winds may potentially introduce discontinuities in the surface wind field that could be manifest in the twin experiment circulation differences as local wind stress curlinduced circulations.Such discontinuities could be eliminated by spatially smoothing the mixed wind fields; however, this was not done here since smoothing may introduce other artifacts into the ocean surface forcing.We instead preferred to preserve the spatial distributions of the wind fields derived from the BHM when present.This aspect of our approach, however, probably deserves further attention.

Data assimilation results
The estimates of Q and Q described in Sect. 4 were used in a sequence of weak constraint 4D-Var data assimilation experiments following the approach outlined in Sect. 2. Two periods were considered corresponding to 2003 and 2005.As noted in Sect.4, twin model experiments and a wind BHM spanning the period 2003-2004 were used to estimate Q and Q which means that weak constraint 4D-Var during any part of this interval will not be independent of the observations since they are used in the twin experiments.Therefore, 2005 represents an independent period, apart from the propagation of information from the data assimilation cycles during 2003 and 2004 into 2005.The 4D-Var configuration is identical to that of N16, only in this case the control vector is augmented with η(t i , t f ), which are the corrections added to the model grid point values to account for model error.The background model error covariance matrix used to compute η is given by W  2. Experiment SNOBHM is the same as NOBHM except that the δx i were spatially smoothed prior to computing to account for the limited sample size.
3. Experiment BHM uses a Q estimated from the twin model runs and the wind BHM as described in Sect.4.3.
4. Experiment SBHM is the same as BHM except that the δ x i were spatially smoothed.
The strong constraint circulation estimates of N16 spanning the same periods were used as a reference, and in all experiments, the strong constraint circulation estimates on 4 January 2003 and 1 January 2005 were used as the background circulation estimates x b on those dates for the first cycle of each sequence of weak constraint 4D-Var calculations.In subsequent data assimilation cycles, the analysis circulation estimate from day 4 of the previous cycle was used as the background estimate for the next cycle using the same 8-day overlapping cycles as described in N16.During each cycle, the corrections η were computed every 6 h, and the realizations of η were assumed to be uncorrelated in time.The η values were interpolated in time to estimate a value at every model time step.While this introduces a correlation in time, the 6 h interpolation interval is very short compared to the 4-day time step dt used to estimate Q.Since the results are quantitatively similar for all four experiments, we will concentrate mainly on experiment BHM and highlight any important differences between experiments where appropriate.
Figure 11 shows time series of the second term on the right-hand side of the cost function J NL in Eq. ( 2) given by J o = (y−H (z a )) T R −1 (y−H (z a )), which is a measure of the weighted difference between the observations and the analysis z a .Since the 4D-Var corrections due to model error are confined to the coast, J o in Fig. 11   rior to the experiment NOBHM, indicating that the BHM is providing useful information about the model error statistics.
The fit of the model to the observations during SBHM and SNOBHM is generally inferior to that of the unsmoothed cases, although still marginally superior to the strong constraint case (not shown).
As discussed in Sect. 4 and Appendix A, strong constraint 4D-Var is likely to compensate for model errors in regions of upwelling and downwelling by making adjustments to the surface forcing fields.When explicit allowance is made for model error during weak constraint 4D-Var, we expect the corrections to the surface forcing to be smaller than those subject to the strong constraint.This is illustrated in Fig. 12d-f, which show the rms differences between the 4D-Var background and analysis estimates of surface wind stress and surface heat flux for the BHM weak constraint 4D-Var calculations during 2005.Also shown in Fig. 12a-c are the corresponding rms differences from the strong constraint 4D-Var calculations of N16. Figure 12a-c indicate that strong constraint 4D-Var makes sizable corrections to the surface stress (∼ 0.02-0.06N m −2 ) and surface heat flux (∼ 150 W m −2 ).The surface freshwater fluxes change little between the strong and weak constraint experiments and are not shown here.In the case of meridional wind stress (Fig. 12b), the 4D-Var corrections are largest near the coast in the vicinity of Cape Mendocino and Cape Blanco.We attribute much of this correction to the influence of model error since COAMPS verifies well against independent observations (Doyle et al., 2009), and as discussed in Sect.3, ROMS CCS is known to possess errors in this region (Broquet et al., 2009a).This is also the same region where the wind BHM indicates that 4D-Var winds are inconsistent with the COAMPS background estimates (cf Fig. 10).Figure 12df reveal that surface forcing corrections during the weak constraint BHM experiment are substantially smaller (∼ 50 %) than those made during the N16 calculations, indicating that 4D-Var places more confidence in the background surface forcing when explicit allowance is made for model errors in the coastal regions.The results are quantitatively similar for the SBHM, NOBHM, and SNOBHM experiments and for The mean weak minus strong differences for 2005 based on experiment BHM are shown in Fig. 13a and c for SST and sea surface salinity (SSS).The mean differences are generally small over much of the model domain and relatively incoherent for SST.In contrast, the standard deviation of the weak minus strong differences for SST and SSS are very coherent as shown in Fig. 13b and d, and reveal that the largest variations occur near the coast, consistent with Fig. 4. Furthermore, the standard deviations are a factor of ∼ 2-6 larger than the mean.The vertical structure of the weak minus strong differences in temperature and salinity is illustrated in Fig. 14 which shows the mean and standard deviation of the differences in each case.In general, the mean differences are relatively incoherent, indicating that the corrections made to the circulation estimates mainly take the form of unbiased errors rather than corrections to account for model bias.In contrast, the patterns of standard deviations are spatially coherent in the upper water column and are generally largest in the vicinity of the thermocline across much of the domain.Therefore, while the corrections for model error are confined mainly to the coastal regions, the impacts on the circulation extend farther offshore.
Recall from Sect. 2 that during weak constraint 4D-Var, corrections are made to account for model error in the time tendency for every prognostic state-vector element at every grid point which correspond to the elements η j of the vector η(t) = {η j (t)}.Figure 15 shows the rms of the corrections η j per time step for the elements j that correspond to SST, SSS, and surface velocity during experiment BHM in 2005.The largest time tendency corrections correspond to the same geographic locations as the largest standard deviations of Q (cf Fig. 4).The vertical structure of η also mirrors that of the standard deviation in Fig. 5, and the corrections η are quantitatively similar for 2003 and for the other experiments (not shown).
The variations of η from one cycle to the next indicate that the model error corrections primarily take the form of unbiased corrections.For example, Fig. 16 shows a time series of the temperature component of η, denoted η T , averaged over each 8-day assimilation cycle, over the entire water column, and over the northern CCS region shown in Fig. 4 during 2003 during experiment BHM.The mean of η T is close to zero (5 × 10 −5 • C per time step), while the standard deviation is significantly larger (1.9 × 10 −4 • C per time step).A nonzero mean would be indicative of corrections of a systematic error.In addition, there is no obvious seasonal dependence in the amplitude of the model error corrections.Time series of other components of η from 2003, 2005, the central CCS, and the other experiments exhibit a similar behavior (not shown).This is a further indication that the weak constraint is correcting for unbiased model errors rather than large systematic errors.

Conclusions
In this paper, we demonstrate how strong constraint data assimilation and a twin model approach can be used in conjunction with linear inverse modeling and a Bayesian hierarchical model to estimate the covariance of model error in a region of coastal upwelling and downwelling.The method has been applied in a state-of-the-art ocean model of the California Current System that is also currently run in near-real time (http://oceanmodeling.ucsc.edu).The CCS is a particularly challenging test region because of the complex nature of the circulation environment which is also characterized by www.adv-stat-clim-meteorol-oceanogr.net/2/171/2016/  a complex system of coastal currents and regions dominated by energetic mesoscale eddies.The estimates of the model error covariance matrix Q were implemented in a weak constraint 4D-Var data assimilation algorithm.Overall, the performance of the weak constraint system is very encouraging in that it indicates that the proposed methodology is able to identify and correct for known deficiencies in the model in the coastal upwelling regions.Furthermore, experiments in which Q is informed by a BHM of surface wind forcing demonstrate the most improvement, in that the fit of the model to the observations is most improved in the coastal regions when compared with strong constraint 4D-Var estimates where the model is assumed to be error-free.Similarly, during the weak constraint experiments, the corrections that are made to the background estimates of surface forcing are greatly reduced compared to those of the strong constraint case, indicating that the data assimilation places more confidence in the surface forcing fields when explicit allowance is made for model errors.4D-Var is based on the assumption that errors in the observations, background, and analysis estimates of the control vector are unbiased and random.The explicit corrections for model error η computed during the weak constraint 4D-Var experiments presented here indicate that these assumptions are generally valid for model error in our experiments.Particularly encouraging is the finding that the Q estimated from a specific time interval (2003)(2004) in the cases considered here) is also effective for correcting for model error during other independent years.
There are, of course, a number of caveats and cautionary notes that should be mentioned.Firstly, the proposed method when applied to surface wind stress increments yields information only about likely model errors in regions where the circulation is tightly coupled to the wind stress, such as regions of upwelling and downwelling.While the signature of model errors on the initial condition increments is potentially more useful for recovering model error information over broader regions, there is currently no practical way to isolate this information.Nonetheless, there are many regions of the world ocean where wind-induced upwelling and downwelling is important, and at the mesoscale the surface wind stress curl and divergence are tightly coupled to the SST.Therefore, the method proposed here could also prove use-  ful in some open ocean regions as well, such as the tropical Pacific.Secondly, the model error covariance estimates computed here are predicated on the signature of model error in the 4D-Var analysis estimates of surface wind stress and the distribution of the background surface wind forcing.However, model errors will also exert undue influence on the surface heat and freshwater fluxes, so an obvious extension of the BHM to compute the posterior distribution of these additional components of the ocean surface forcing is also warranted.Thirdly, the model error estimates derived here will depend on the prior constraints imposed on the data assimilation system via the background error covariance B, the observation error covariance R, and the distribution of the observations.For example, in areas that are devoid of observations, the surface forcing corrections made by 4D-Var may be small.This does not necessarily mean that model error is absent in such regions, only that there is no direct information about the presence (or lack thereof) of such errors.However, further analysis (not shown) indicates that the rms surface flux differences between the 4D-Var analyses and background fields in Fig. 12 are not obviously related in any way to the distribution of the observations during the same period.Therefore, we feel confident that the 4D-Var surface forcing corrections are not being overly influenced by the observation sampling in the present experiments.Finally, it is useful to speculate on the factors that may be contributing to model errors in coastal regions in the present study.A very obvious omission in the ROMS CCS configuration used here are the sources of freshwater associated with the Columbia River and the Juan de Fuca Strait.This may account in part for the geographical distribution of the largest model error covariances along the coast of Oregon and Washington (cf.Fig. 4) and the freshening of coastal waters during weak constraint 4D-Var (cf.Fig. 13c).One could argue that this represents a source of forcing error rather than model error.However, as the results of this study indicate, the nature of the corrections for model error are more consistent with that of an unbiased error rather than of a correction for a bias such as one might expect in the presence of a persistent error in the freshwater flux.Another conspicuous omission in the current configura-tion of the model is tidal forcing.Intensification of tidal circulations can occur in the vicinity of coastal bathymetry and topography along the US west coast (e.g., Osbourne et al., 2014), and may represent another potential source of error consistent with the localized nature of Q in Fig. 4. Other potential sources of model error include bathymetry.The continental shelf along the US west coast is generally narrow and not well resolved by the 10 km model grid used here.In addition, the required smoothing of the bathymetry in terraincoordinate-following models introduces further errors.There are several bathymetric features, such as Heceta Bank, that influence the circulation in the northern reaches of the CCS which are poorly resolved in the present model and consistent with the geographic structure of the variance of Q in Fig. 4.

Data availability
The total volume of data processed in the experiments reported here exceeds 15 TB, so it is not feasible to make them available online.The strong constraint analysis of N16 is, however, available at http://oceanmodeling.ucsc.edu/reanalccs13.Subsets of the data may be requested by contacting the corresponding author directly.
Under the assumption of the strong constraint, t (t) will have an expression on x(t i ), f τ and f Q during S1 while during S2 there will be an additional contribution from r (t).The difference between corresponding members of S1 and S2 will be denoted x and will be associated primarily with the influence of only r (t) on the elements of the control vector.Therefore, by subtracting Eq. (A1) from Eq. (A2), the influence of the actual model error t (t) can be temporarily eliminated (assuming t (t) and r (t) are independent), and to first order the differences x = x S2 − x S1 will be described by ) where N| x S1 denotes the tangent linear model linearized about the analyses from sequence S1.Since the prior initial conditions are identical for corresponding cycles of each 4D-Var sequence, the tangent linear Eq. (A3) can be used to cleanly separate the influence of each component of the control vector (i.e., x(t i ), f τ and f Q ) on the circulation estimates.For example, if we denote x i as the difference between the analyses of S2 and S1 that are associated with differences in x(t i ) alone, we can use Eq.(A3) to explore the time evolution of these differences.Similarly, we can use Eq.(A3) to quantify the analysis differences x τ due to differences in f τ alone and x Q h due to f Q alone.The time series of x averaged over each 8-day 4D-Var cycle were then modeled as a first-order Markov process as described in Sect.4.2 (cf Eq. 12) to determine the degree to which the covariance properties of r (t) can be recovered from the differences in different elements of the control vector arising from strong constraint 4D-Var.
The experiment S2 was performed using time series of r (t) drawn from the distribution N(0, Q).As noted in Sect.4.2, the background model error covariance matrix can be expressed as , where is the matrix of EOFs of the δx and Q is the covariance matrix of the stochastic forcing of the associated principal components.Furthermore, Q can be written as E E T , where E is the matrix of eigenvectors of Q and is the diagonal matrix of eigenvalues.The multivariate random vector r (t) = E 1 2 ξ (t) represents a random draw from a Gaussian probability distribution with covariance Q, where ξ ∼ N(0, I).Two different representations of Q were considered here.In the first case, Q = Q 1 was based on the identified in Sect.4.3 and describes errors primarily in the coastal upwelling and downwelling regions.In the second case, Q = Q 2 was constructed from a random sample of 50 spatially smoothed analysis increments from a randomly chosen year (2008) of N16.From the sample of 50 increments, a new orthonormal basis was constructed using a Gram-Schmidt procedure.For convenience, the eigenspectrum of Q 2 was chosen to be proportional to that of Q 1 and rescaled to give a similar total variance for the two estimates of r (t).
Figure A1 shows a summary of the results of the experiments using Q 1 in terms of the standard deviation σ sst of the artificial "model error" in SST added to ROMS.These results are representative of other fields also.The expected error associated with Q 1 for SST is shown in Fig. A1a and highlights the coastally trapped nature of r (t). Figure A1b shows an estimate of σ sst derived from the first-order Markov model of x i , the circulation associated with the initial condition differences resulting from r (t).Clearly the additional artificial "model error" can be quite successfully recovered from the signature of r in x(t i ), even given the relatively small sample size used here (92 realizations of x). Figure A1c shows σ sst derived from x τ , the circulation associated with the surface wind stress differences resulting from r (t).Clearly, the wind stress is also able to recover the pattern of σ sst quite well, although the amplitude is somewhat lower than the expected value.Although when x only on day 8 is used, the estimates of σ sst are significantly larger (not shown).The surface heat and freshwater fluxes contribute little to the estimates of σ sst and are not shown.Figure A1a-c demonstrate the important point that there is a unique rela- tionship between the expression of the model error r in the upwelling circulation and the expression of model error in the surface wind stress increments arising from strong constraint 4D-Var, the basic hypothesis on which the method presented in Sect. 4 is based.
Figure A1d-f presents a summary of σ sst for the experiments using Q 2 .Figure A1d shows the expected artificial model error variance for SST, which, in contrast to the previous case, spans the entire model domain.Figure A1e shows σ sst derived from x i , and as in the case of Q 1 , the total variance in SST due to model error can be recovered fairly well from the first-order Markov model for x i alone.Con-versely, Fig. A1f derived from x τ shows that r in the deep ocean has little expression in f τ .Only in the regions of coastal upwelling and downwelling does r have any appreciable expression on f τ , and the first-order Markov models of x τ can only recover the coastal component of Q 2 .The same is true for x Q associated f Q (not shown).
Figure A1 suggests that the model error covariance can be reliably estimated everywhere from the expression of model error in all elements of the control vector in combination, but particularly from the initial conditions.However, an operational practical environment is quite different from the nature of the experiments presented here, where we chose the prior www.adv-stat-clim-meteorol-oceanogr.net/2/171/2016/ initial conditions to be invariant for a given cycle start date.In an operational environment, we will never be able to identify the component of x(t i ) that is associated solely with the model error.This would require an estimate of the probability distribution of x(t i ) based on which we could decide which corrections to x b (t i ) are legitimate corrections for errors in the background initial conditions and which corrections are most likely the result of model errors.For the wind stress f τ , however, an estimate of the probability distribution in the form of a BHM is available as described in Sect.4.3.Construction of a BHM for the surface wind stress is tractable because of the relatively low dimension of the problem and the availability of good observation coverage from scatterometers.The same approach would be intractable for x(t i ), however, because of the multivariate nature of the state vector and the much larger dimension of the problem and because of the lack of high density observations for the entire ocean state.If the latter were available, we would probably not need ocean data assimilation at all.
The experiments presented here explore how the signature of model error is manifested in the surface wind stress over an 8-day assimilation cycle.It seems reasonable to expect that the influence of model error on the surface winds will increase with the length of the assimilation window, so there may be value in considering windows longer than 8 days in order to better quantify model error.However, care must be exercised using this approach since the tangent linear assumption on which 4D-Var is predicated will be violated if the assimilation window becomes too long.
, ROMS has three open boundaries that were constrained by time-evolving circulation fields from the global Simple Ocean Data Assimilation (SODA; version SODA POP 2.2.4) product of Carton and Giese (2008).These constitute the background lateral boundary conditions b b .

Figure 1 .
Figure 1.The model domain and bathymetry used in the present study.A schematic representation of some of the important circulation features is also shown.

Figure 2 .
Figure 2. (a) An estimate of the leading 120 eigenvalues λ of C 0 based on 120 iterations of the Lanczos algorithm using unsmoothed (red curve) and smoothed (black) δx i .(b) Log 10 of formal error estimates of | C 0 ψ j − λψ j | vs. EOF number j using unsmoothed (red curve) and smoothed (black) δx i .(c) The cumulative percent variance explained vs. the number of EOFs using unsmoothed (red curve) and smoothed (black) δx i .

Figure 3 .
Figure 3.The autocorrelation vs. lag for the stochastic forcing associated with each principle component element of p in Eq. (14).

Figure 4 .
Figure 4.The standard deviation derived from Q (computed using unsmoothed δx) associated with the surface fields of (a) temperature (C per dt), (b) salinity (per dt), (c) current speed (m s −1 per dt), and (d) sea surface height (m per dt).The black lines delineate a central and northern CCS region used to compute several diagnostics in Sect. 5.

Figure 6 .
Figure 6.(a) The eigenspectrum of Q derived from smooth (red) and unsmoothed (black) δx.(b) The cumulative percentage variances explained by the EOFs of Q.

Figure 7 .
Figure 7.The structure of the leading EOF of Q (computed using unsmoothed δx) for surface (a) temperature, (b) salinity, (c) meridional velocity, and (d) sea surface height.

Figure 9 .
Figure 9. CCS surface wind BHM posterior distribution summaries for 19 January 2003.Velocity component variances (m 2 s −2 ) for (a) u and (b) v as estimated in the BHM posterior distribution.Uncertainties, expressed as variances, arise from data stage coverage (swath patterns are evident in a and b) and from process model misfits.The spread in velocity vector clusters (red) in (c) provides an intuitive sense of these uncertainties as well.In (c), 10 surface wind realizations (red) have been randomly selected from the posterior distribution (e.g., see Milliff et al., 2011) in the same subdomain for Fig. 8. Posterior wind vectors from the 4D-Var are shown in black in panel (c).The strongconstraint 4D-Var corrections have incrementally adjusted the surface wind forcing to unrealistic values outside the distribution obtained from the BHM for this region (see text).

Figure 10 .
Figure 10.The root mean square (rms) difference between the COAMPS background surface wind stress and the "mixed wind" stress for the period 2003-2004.The rms differences are shown as vectors for more clarity.
4. Four sequences of weak constraint 4D-Var circulation estimates were computed for 2003 and 2005: 1. Experiment NOBHM uses a Q estimated from the twin model runs alone, making no use of the posterior wind distributions from the BHM.

Figure 11 .
Figure 11.Time series of the observational component of the cost function, J o , for observations within the two coastal regions shown in Fig. 4 that delineate a region spanning the northern CCS and the central CCS to 100 km offshore during 2003 (a and b) and 2005 (c and d).The regional J o is shown for the strong constraint case (blue curve), and the weak constraint experiments BHM (red curve) and NOBHM (green curve).
was computed only for observations that fall within the two regions adjacent to the coast shown in Fig. 4 that span the northern and central CCS and extend 100 km offshore.Time series of these regional J o values are shown for the strong constraint case of N16 and the weak constraint experiments BHM and NOBHM for both 2003 and 2005.The posterior fit of the model to the observations is similar in all cases, but is generally lowest during experiment BHM indicating that accounting for model error brings the circulation estimates closer to the observations.Much of this improvement is associated with a better fit of the model to the satellite observations (not shown).The change in the degree of fit of the model to the observations is not expected to change dramatically between the strong and weak constraint experiments since the strong constraint system is already performing well (see N16).Nonetheless, Fig. 11 also shows that experiment BHM is generally supewww.adv-stat-clim-meteorol-oceanogr.net/2/171/2016/Adv.Stat.Clim.Meteorol.Oceanogr., 2, 171-192, 2016

Figure 12 .
Figure 12.The rms difference between analysis and background 4D-Var estimates of (a) zonal wind stress (N m −2 ), (b) meridional wind stress, (N m −2 ), and (c) surface heat flux (W m −2 ) from the strong constraint calculations of N16 for 2005.Panels (d)-(f) show the corresponding differences for the weak constraint 4D-Var experiment BHM.

Figure 13 .
Figure 13.The mean weak minus strong constraint 4D-Var differences for (a) SST and (c) SSS for the experiment BHM estimates during 2005.The standard deviation for the (b) SST and (d) SSS weak minus constraint differences are also shown.

Figure 14 .
Figure 14.Vertical sections of the mean weak minus strong constraint 4D-Var differences for (a) temperature and (c) salinity for the experiment BHM estimates during 2005.The standard deviation for the (b) temperature and (d) salinity weak minus strong constraint differences are also shown.

Figure 15 .
Figure 15.The rms corrections η made by weak constraint 4D-Var per time step during experiment BHM in 2005 for (a) SST, (b) SSS, (c) surface zonal velocity, and (d) surface meridional velocity.

Figure 16 .
Figure 16.Time series of the cycle-average temperature components of η averaged over the entire water column and over the northern CCS region indicated in Fig. 4 during 2003 for experiment BHM.The units are • C per time step.

Figure A1 .
Figure A1.The standard deviation of the additional artificial source of model error added to experiments S2: (a) the expected error based on Q 1 , (b) derived from x i using Q 1 during S2, and (c) derived from x τ using Q 1 during S2.(d) The expected error based on Q 2 , (e) same as (b) but using Q 2 during S2, and (f) same as (c) but using Q 2 during S2.The color bar on the left corresponds to panels (a), (b), (d), and (e), while the color bar on the right is for (c) and (f).

Recovering model errors from strong constraint 4D-Var forcing corrections
The underlying premise of the method used in Sect. 4 to recover an estimate of Q is that estimates of additive model errors can be obtained from the influence that such errors have on ocean surface forcing fields in regions of coastal upwelling and downwelling and for which estimates of the probability distribution are available.The efficacy and limitations of this premise are demonstrated here using two sets of 4D-Var experiments in which artificial model errors were added to ROMS.Following this approach, two sequences of strong constraint 4D-Var experiments for 2003 using 8day cycles were performed comprising three sets of runs for each using the same configuration as N16.During the first sequence, denoted S1, 4D-Var was run adjusting the initial conditions (x(t i )), surface wind stress (f τ ), and surface heat and freshwater fluxes (f Q ).The second sequence, denoted S2, was the same as S1 but included an additional artificial source of model error r (t) with known covariance imposed on the model.The start days t i for each cycle were the same in both experiments.The prior initial conditions x b (t i ) varied with t i but were the same in both experiments for a given t i .During S1, the 4D-Var analyses are influenced by the actual model t (t), which is unknown, while during S2 they are influenced by both t (t) and r (t).