# Some Aspects of Sensitivity Analysis in Variational Data Assimilation for Coupled Dynamical Systems

Authors: Sergei Soldatenko; Peter Steinle; Chris Tingwell; Denis Chichkine
Journal: Advances in Meteorology
 (2015)

---

## Abstract

Variational data assimilation (VDA) remains one of the key issues arising in many fields of geosciences including the numerical weather prediction. While the theory of VDA is well established, there are a number of issues with practical implementation that require additional consideration and study. However, the exploration of VDA requires considerable computational resources. For simple enough low-order models, the computational cost is minor and therefore models of this class are used as simple test instruments to emulate more complex systems. In this paper, the sensitivity with respect to variations in the parameters of one of the main components of VDA, the nonlinear forecasting model, is considered. For chaotic atmospheric dynamics, conventional methods of sensitivity analysis provide uninformative results since the envelopes of sensitivity functions grow with time and sensitivity functions themselves demonstrate the oscillating behaviour. The use of sensitivity analysis method, developed on the basis of the theory of shadowing pseudoorbits in dynamical systems, allows us to calculate sensitivity functions correctly. Sensitivity estimates for a simple coupled dynamical system are calculated and presented in the paper. To estimate the influence of model parameter uncertainties on the forecast, the relative error in the energy norm is applied.


---

## Body

## 1. Introduction

The earth system consists of several interactive dynamical subsystems and each of them covers a broad temporal and spatial spectrum of motions and physical processes. The components of the earth system have many differences in their physical properties, structure, and behavior but are linked together by fluxes of momentum and mass as well as sensible and latent heat. All of these subsystems interact with each other in different ways and can be strongly or weakly coupled. Prediction of future state of the earth system and its components is one of the most important problems of modern science. The most significant progress has been achieved in the forecasting of the atmosphere via numerical models, which describe the dynamical and physical processes in the earth’s gaseous envelope. It is clear that further improvement of forecasts can be pursued via the development of coupled modeling systems that primarily combine the atmosphere and the ocean and describe the interactions between these two systems. Since numerical weather prediction systems calculate a future state of the atmosphere and ocean by integrating a set of partial differential equations that describe the fluid dynamics and thermodynamics, initial conditions that accurately represent the state of the atmosphere and ocean at a certain initial time must be formulated. Numerical weather prediction systems use data assimilation procedures to estimate initial conditions for forecasting models from observations. Data assimilation remains one of the key issues not only in the numerical weather prediction (NWP) but also in other geophysical sciences.One of the most advanced and effective data assimilation techniques is four-dimensional variational data assimilation (4D-Var). In particular, the weather forecasts produced by the ACCESS (Australian Community Climate and Earth System Simulator) at the Bureau of Meteorology use 4D-Var in the incremental formulation developed at the Met Office [1, 2]. A comprehensive historical review and current status of 4D-Var are presented, for example, in [3–7]. In our view point, some papers and books, such as [8–17], have contributed significantly to the development of mathematical foundation, theory, and practice of 4D-Var.In general, the main objective of 4D-Var is to define, as perfectly as possible, the state of a dynamical system by combining, in statistically optimal manner, the observations of state variables of a real physical system together with certain prior information. In NWP this prior information is usually referred to as the background. Mathematically, 4D-Var procedures are formulated as an optimization problem, in which the initial condition plays the role of control vector and model equations are considered as constraints. While the theory of variational data assimilation is well established, there are a number of issues with practical implementation that require additional consideration and study. The performance of 4D-Var schemes depends on their key information components, such as the available observations, estimates of the observation, and background error covariances that are quantified by the corresponding matrices, as well as the background state. All of those information components strongly impact the accuracy of calculated initial conditions, thus influencing the forecast quality. Thus, it is important to estimate the sensitivity of a certain forecast aspect with respect to variations in the observational data, background information, and error statistics. This problem is usually formulated within the framework of the adjoint-based sensitivity analysis (e.g., [18–29]).However, an adjoint model used to calculate sensitivity functions is derived from a linearized forward propagation model and contains numerous input parameters. Consequently, linearization of strongly nonlinear NWP models and also uncertainties in their numerous parameters generate errors in the initial conditions obtained by data assimilation systems. The influence of linearization and parameter uncertainties on the results of data assimilation can be studied, ideally for each particular NWP model, using sensitivity analysis [30, 31]. Even more problems arise when considering coupled 4D-Var data assimilation schemes since the atmosphere and ocean have very different physical properties and time-space spectrum of motions generating initialization shock. Several coupling strategies are being developed for use in NWP systems; however, all of them introduce issues that require additional detailed consideration. For example, the influence of coupling strength on the initial conditions obtained by 4D-VAR procedures is one such issue that is important to study and analyze.Good practice in the development of NWP models and data assimilation systems requires evaluating the confidence in the model. In this context, it is important to estimate the influence of parameter variations on system dynamics and to find those parameters that have the largest impact on system behaviour. Sensitivity analysis, which is an essential element of model building and quality assurance, addresses this very important issue.The exploration of coupled 4D-Var systems, parameter estimation, and sensitivity analysis require considerable computational resources. For simple enough low-order coupled models, the computational cost is minor and, for that reason, models of this class are widely used as simple test instruments to emulate more complex systems. In this paper, we describe a coupled nonlinear dynamical system, which is composed of fast (the “atmosphere”) and slow (the “ocean”) versions of the well-known Lorenz [32] model. This low-order coupled system allows us to mimic the atmosphere-ocean system and therefore serves as a key element of a theoretical and computational framework for the study of various aspects of coupled 4D-Var procedures [33, 34]. Under certain conditions the Lorenz model exhibits a chaotic behaviour and using conventional methods of sensitivity analysis can be questionable in terms of interpretation of the obtained results [35–37]. The “shadowing” method [36, 37] for estimating the system sensitivity to variations in its parameters allows us to calculate the average along the trajectory sensitivities and therefore to make a clear conclusion with respect to the system sensitivity to its parameters. This method is based on the pseudoorbit shadowing in dynamical systems [38, 39]. Calculated sensitivity coefficients obtained via conventional methods and the “shadowing” approach are presented in the paper.We also succinctly consider commonly used techniques for sensitivity analysis and parameter estimations of dynamical systems and study the influence of coupling strength parameter on the dynamical behaviour of the coupled system using Lyapunov characteristic exponent analysis. It was found that the coupling strength parameter strongly affects the system dynamics both quantitatively and qualitatively. This fact should be taken into consideration when choosing the coupling strategy in data assimilation systems.


## 2. Low-Order Coupled Dynamical System

In this section we consider a low-order coupled nonlinear dynamical system obtained by coupling of two versions of the original Lorenz model (L63) [32] with distinct time scales, which differ by a factor ε (e.g., [33, 34]):(1a)x˙=σy-x-caX+k,y˙=rx-y-xz+caY+k,z˙=xy-bz+czZ,(1b)X˙=εσY-X-cx+k,Y˙=εrX-Y-aXZ+cy+k,Z˙=εaXY-bZ-czz,where lower case letters represent the fast subsystem and capital letters the slow subsystem, σ, r, and b are the parameters of L63 model, c is a coupling strength parameter for the x and y variables, cz is a coupling strength parameter for z, k is an “uncentering” parameter, and a is a parameter representing the amplitude scale factor. The value a=1 indicates that two systems have the same amplitude scale. Thus, the state vector of the coupled model (1a) and (1b) is x=x,y,z,X,Y,ZT and the model parameter vector is α=σ,r,b,a,c,cz,k,εT. Without loss of generality, we can assume that a=1, k=0, and c=cz, then α=σ,r,b,c,εT and the system (1a) and (1b) can be rewritten in the operator form:(2)dxdt=L+Qx,where the nonlinear uncoupled operator L and linear coupled operator Q are represented by the following matrices:(3)L=-σσ0r-1-x00x-b-εσεσ00εr-ε-εX0εX-εb,Q=-c0000c000c-c000c0000-c.The unperturbed parameter values are taken as(4)σ0=10,r0=28,b0=83,ε0=0.1,c0∈0.1;1.2.Chosen values of σ,r, and b correspond to chaotic behaviour of the L63 model. For σ=10 and b=8/3 the critical value of parameter r is 24.74, which means that any value of r larger than 24.74 induces chaotic behaviour [32]. The parameter ε=0.1 indicates that the slow system is 10 times slower than the fast system.The coupling strength parameterc plays a very important role in qualitative changes in the system dynamics since this parameter controls the interactions between fast and slow subsystems. Qualitative changes in the dynamical properties of a system can be detected by determining and analyzing the system’s spectrum of Lyapunov exponents, which characterize the average rate of exponential divergence (or convergence) of nearby trajectories in the phase space. In the analysis of coupled dynamical systems we are dealing with conditional Lyapunov exponents that are normally used to characterize the synchronization with coupled systems. System (2) has six distinct exponents. If the parameter c tends to zero, then system (2) has two positive, two zero, and two negative Lyapunov exponents. The influence of coupling strength parameter c on the two largest conditional Lyapunov exponents is illustrated in Figure 1. The numerical experiments demonstrated that those, initially positive, exponents decrease monotonically with an increase in the parameterc. At about c≈0.8 they approach the x-axis and at about c≈0.95 negative values. Thus, for c>0.95, the dynamics of both fast and slow subsystems become phase-synchronous [40]. For c>1.0, a limit circle dynamical regime is observed since all six exponents become negative.Figure 1
Two largest conditional Lyapunov exponents as functions of coupling strength parameter.Apart from Lyapunov exponents, autocorrelation functions (ACFs) enable one to distinguish between regular and chaotic processes and to detect transition from order to chaos. In particular, for chaotic motions, ACF decreases in time, in many cases exponentially, while for regular motions, ACF is unchanged or oscillating. In general, however, the behaviour of ACFs of chaotic oscillations is frequently very complicated and depends on many factors (e.g., [41, 42]). Autocorrelation functions can also be used to define the so-called typical time memory (typical timescale) of a process [43]. If it is positive, ACF is considered to have some degree of persistence: a tendency for a system to remain in the same state from one moment in time to the next. The ACF for a given discrete dynamic variable xii=0N-1 is defined as(5)Cs=xixi+s-xixi+s,where the angular brackets denote ensemble averaging. Assuming time series originates from a stationary and ergodic process, ensemble averaging can be replaced by time averaging over a single normal realization:(6)Cs=umum+s-u2.Signal analysis commonly uses the normalized ACF, defined as Rs=Cs/C0. ACF plots for realizations of dynamic variables x and X, and z and Z calculated for different values of the coupling strength parameter c are presented in Figures 2 and 3, respectively. For relatively small parameter c   (c<0.4), the ACFs for both x and X variables decrease fairly rapidly to zero, consistently with the chaotic behaviour of the coupled system. However, as expected, the rate of decay of the ACF of the slow variable X is less than that of the fast variable x. The ACF envelopes for variables z and Z also decay almost exponentially from the maximum to zero. For coupling strength parameter on the interval 0.4<c<0.6 the ACF of the fast variable x becomes smooth and converges to zero. At the same time, the envelopes of the ACFs of variables X, z, and Z demonstrate a fairly rapid fall, indicating the chaotic behaviour. As the parameter c increases, the ACFs become periodic and their envelopes decay slowly with time, indicating transition to regularity. For c>0.8 calculated ACFs show periodic signal components.Figure 2
Autocorrelation functions for dynamic variablesx and X for different parameter c.Figure 3
Autocorrelation functions for dynamic variablesz and Z for different parameter c.


## 3. The Basics of Four-Dimensional Variational Data Assimilation

Atmospheric models used for operational NWP are mainly deterministic and derived from a set of multidimensional nonlinear differential equations in partial derivatives, which are the equations of fluid dynamics and thermodynamics that describe atmospheric processes and atmosphere-underlying surface interactions. A generic atmospheric model can be represented by the following continuous autonomous dynamical system:(7)dxtdt=fxt,αt,t∈0,τ=T,x0=x0.Here xt is the state vector belonging to a Hilbert space X, α∈P is a parameter vector, where P is a Hilbert space (the space of the model parameters), f is a nonlinear vector-valued function, such that f:X×P×T→X, andx0 is a given vector function. This infinite-dimensional model has to be truncated by some means to finite-dimensional approximate model, for which a solution can be sought numerically. Applying either a projection onto a finite set of basic functions or a discretization in time and space, one can derive the discrete atmospheric model which can be represented as s discrete nonlinear dynamical system given by the equation(8)xi+1=Mi,i+1xi+εim,where xi∈Rn is the n-dimensional state vector representing the complete set of the model variables that determine the internal state of the atmospheric model at time ti, Mi,i+1 is the nonlinear operator that indirectly contains model parameters and propagates the state vector from time ti to time ti+1 for i=0,…,N-1, and εim is the model errors. It is usually assumed that model (8) is “perfect” (εim=0), that is the forecast has no errors if the initial condition is perfect. In this case, given the model operator and the initial condition x0, (8) uniquely specifies the orbit of the dynamical system. Let yi0∈Rm be the m-dimensional vector of observations measured at a discrete time ti, i=0,…,N that are linked to the system state via the following equation:(9)yi0=Hixi+εi0,where Hi:Rn→Rm is the nonlinear observation operator that maps the state vector to observation space. It is usually assumed that the observation errors εi0 are unbiased, serially uncorrelated, and normally distributed with known covariance matrices Ri∈Rm×m.Suppose that at the initial time the prior (background) model statex0b is known and represents the “best” estimate of the “true” state x0t before any observations are taken. This background state is provided by a previous forecast. It is assumed that x0b has unbiased and normally distributed errors εb with known covariance matrix B0∈Rn×n:(10)x0b=x0t+εb.Given the observationsyio at time ti, the corresponding observation error covariance matrices Ri (i=0,…,N), the background initial state x0b, and the error covariance matrix B0, the 4D-Var data assimilation seeks to minimize, with respect to x0, a certain cost function Jx expressing the “distance” between observations and corresponding model state using the model equations as constraints:(11)x0a=arg minJxsubject to x satisfying the set of the “ideal” (εim=0) model equations with initial state x0:(12)xi+1=M0,i+1x0,i=1,…,N.The 4D-Var cost function is usually written as (e.g., [17, 18]):(13)Jx0=12x0-x0bTB0-1x0-x0b+12∑i=0NHixi-yi0TRi-1Hixi-yi0.The optimization problem (11) is nonlinear with strong constraints and an iterative minimization algorithm (e.g., gradient-based technique) is required to obtain the solution. The gradient of the cost function (13) is as follows:(14)∇x0Jx0=B0-1x0-x0b+∑i=0NM0,iTHiTRi-1Hi(xi)-yi0,where M0,iT is the adjoint of the linearized model operator M0,i=M0,i′xi and HiT is the adjoint of the linearized observation operator Hi=Hi′xi. If the model is not “perfect” then we need to take into account the model errors εm, which are sometimes taken as Gaussian noise:(15)εm∈N0,Q,where Q is a model error covariance matrix. Thus, we obtain the weakly constrained 4D-Var data assimilation and the following term(16)Jm=12∑i=1Nxi-Mi-1,ixi-1TQi-1xi-Mi-1,ixi-1should be added to the right-hand side of the cost function (13) [6, 19].It is important to make the following comments. While NWP has significantly improved over the last several decades, weather forecasts are still intrinsically uncertain. This is because mathematical models used in NWP have several sources of uncertainty and therefore a number of sources of errors. The first one is an intrinsic uncertainty due to chaotic nature of the system. The second one is a structural uncertainty, which is how the model itself represents the physical processes incorporated into it. It is important to underline that our knowledge of the earth system is always imperfect and, therefore, we can only theoretically design the “ideal” model. However, improving model physics demonstrated that even theoretically well-posed models failed to accurately simulate and predict the dynamics of real atmospheric processes. This is because numerical models have a parametric uncertainty (how accurate model parameters are) initial and boundary uncertainty (are initial and boundary conditions known precisely) and, in addition, numerical errors. All of those uncertainties (errors) impair the weather forecast accuracy and limit the time horizon of accurate NWP. Current time horizon of synoptic-scale NWP is several days.The climate study has significantly longer time horizon: several decades. The sensitivity analysis of the climate system is associated with stability of characteristics of climatic model attractors with respect to perturbations in model parameters.A priori estimation of the behaviour of state vector, when the perturbations in the climate model parameters tend to zero, is generally an unresolved problem since it is not known whether or not the invariant measure of climate modeling system is continuous with respect to the small perturbations in the differential matrix operator of the numerical model. Indeed, at certain model parameter values, different bifurcation can occur in the system phase space. Therefore, dynamics on the attractor generated by the model may change considerably even for small parameter perturbations.


## 4. Sensitivity Analysis: Essentials of Conventional Approaches

One of the commonly used measures for estimating the influence of model parameter variations on the state variables is the sensitivity coefficient, which is the derivative of a certain component of a model state vectorxi with respect to some model parameter αj [30, 31, 44]:(17)Sij≡∂xi∂αjαj=αj0=limδα→0xiαj0+δαj-xiαj0δαj,where δαj is the infinitesimal perturbation of parameter αj around some fixed point αj0. Differentiating (8) with respect to α, we obtain the set of nonhomogeneous ODEs, the so-called sensitivity equations, which can be written as(18)dSjdt=M·Sj+Djj=1,…,m,where Sj=∂x/∂αj=S1,j,S2,j,…,Sn,jT is the sensitivity vector with respect to parameter αj, Dj=∂f1/∂αj,∂f2/∂αj,…,∂fn/∂αjT, and M is a Jacobian matrix. Thus, to analyze the sensitivity of system (7) with respect to parameter αj one can solve the following set of differential equations with given initial conditions:(19)dxdt=fx,α,x0=x0,dSjdt=M·Sj+Dj,Sj0=Sj0.Sensitivity equations describe the evolution of sensitivity coefficients along a given trajectory and therefore allow tracing the sensitivity dynamics in time. The procedure for computing sensitivity coefficients includes the following steps.(1)
Obtain initial conditions on the system attractor at timet=t0 by integrating the nonlinear model equations (7) for a long enough time range t0,τ, starting from random initial conditions.(2)
Solve the nonlinear model equations (7) to calculate a trajectory xt, t∈t0,τ.(3)
Calculate a model Jacobian matrix and a parametric Jacobian matrix.(4)
Solve the sensitivity equations (18) with given initial conditions to obtain the desired sensitivity coefficients.Sensitivity analysis allows us also to explore the sensitivity of a generic objective function (performance measure), which characterizes the dynamical system (7):(20)Jx,α=∫0τΦt;x,αdt,where Φ is a nonlinear function of the state variables x and model parameters α. The gradient of the functional J with respect to the parameters α around the unperturbed state vector x0(21)∇αJx0,α0=dJdα1,…,dJdαmTx0,α0quantifies the influence of parameters on the model output results. In particular, the effect of the jth parameter can be estimated as follows:(22)dJdαjαj0≈Jx0+δx0;α10,…,αj0+δαj,…,αm0-Jx0,α0δαj,where δαi is the variation in parameter αi0. Note that(23)dJdαj=∑i=1n∂J∂xi∂xi∂αj+∂J∂αj=∑i=1nSij∂J∂xi+∂J∂αj.This approach is acceptable for low-order models. However, the accuracy of sensitivity estimates strongly depends on choice of the perturbation δαi. By introducing the Gâteaux differential, the sensitivity analysis problem can be considered in the differential formulation eliminating the need to set the value of δαi [30, 31]. The Gâteaux differential for the objective function (20) has the following form:(24)δJx0,α0;δx,δα=∫0τ∂Φ∂xx0,α0·δx+∂Φ∂αx0,α0·δαdt.Here δx is the state vector variation due to the variation in the parameter vector in the direction δα. Linearizing the nonlinear model (7) around an unperturbed trajectory x0t, we obtain the following system of variational equations, the so-called tangent linear model, for calculating δx:(25)∂δx∂t=∂f∂xx0,α0·δx+∂f∂αx0,α0·δα,t∈0,τ,δx0=δx0.Then using (24) we can calculate the variation δJ. Since δJx0,α0;δx,δα=∇αJ,δα, where ·,· is a scalar product, the model sensitivity with respect to parameter variations can be estimated by calculating the components of the gradient ∇αJ. However, this method is computationally ineffective if the number of model parameters m is large. The use of adjoint equations allows obtaining the required sensitivity estimates within a single computational experiment (e.g., [6, 30, 31]) since the gradient ∇αJ can be calculated by the following equation:(26)∇αJx0,α0=∫0τ∂Φ∂αx0,α0-∂f∂αx0,α0T·x∗dt,where the vector function x∗ is the solution of adjoint model(27)-∂x∗∂t-∂f∂xx0,α0Tx∗=-∂Φ∂xx0,α0,t∈0,τ,x∗τ=0.This equation is numerically integrated in the inverse time direction. Thus, the algorithm for computing sensitivity functions is as follows.(1)
Obtain initial conditions on the system attractor at timet0 by integrating the nonlinear model equations (7) for a long enough time range t0,τ, starting from random initial conditions.(2)
Solve the nonlinear model equations (7) to calculate a trajectory xt, t∈t0,τ.(3)
Calculate the right-hand side of (27) and then integrate numerically this equation in the inverse time direction with the initial conditions x∗τ=0.(4)
Calculate the gradient (26).


## 5. Testing the Variational Data Assimilation

We will consider the dynamics of system (2) on its attractor. In order to obtain the model attractor, the numerical integration of (2) is started at tD=-20 with the initial conditions (28)x(tD)=0.01,0.01,0.01,0.02,0.02,0.02Tand finished at t0=0 to guarantee that the calculated model state vector x0=x(0) is on the model attractor. The forecast, obtained by NWP models, is substantially determined by initial conditions, which calculated via 4D-Var. The accuracy of initial conditions strongly depends on numerous parameters of 4D-Var systems. Some of these parameters are uncertain. To estimate the influence of parameter variations on the forecast obtained by model (2), the “true” and “forecast” trajectories were calculated. The “true” trajectory xtti and the “true” state xttf at the verification time tf are obtained by integrating the model equations over the time interval t0,tf with unperturbed parameter vector α0 and initial conditions x0. Then, the forecast trajectory xftk and the forecast state xftf at tf are obtained by integration of the “forecast” model (2) with initial conditions x0 and a certain perturbed model parameter. Thus, the “forecast” model has the same set of equations as the “true” model; however, some model parameter is slightly changed (i.e., this parameter is known with uncertainty). In order to measure forecast errors the relative error in energy norm is used: (29)exf=xf-xtTxf-xtxtTxt1/2.As an example, variations in parametersr and c are considered and forecast is made for two time periods tf: 2.5 and 5.0 of nondimensional time units. The forecast error (29) is estimated at time tf. Table 1 shows the results of forecast verifications. It is obvious that the less the forecast error measure exf, the higher the forecast skill. Qualitatively, calculated results are consistent with real numerical weather forecasts obtained with complex state-of-the-art NWP models: longer tf leads to lower forecast accuracy and smaller parameter variations (difference between the “true” parameter value and the value used in the “real” model) lead to better forecast accuracy. It can also be observed that parameter r influences the forecast accuracy almost twice as much as parameter c.Table 1
Relative errors in energy norm (subscript at parameters denotes the nondimensional time unit).
Parameter variations in percent of unperturbed value
1%
−1%
5%
−5%
10%
−10%
r
2.5
0.1273
0.1266
0.2871
0.2796
0.4081
0.3890
c
2.5
0.0584
0.0583
0.1312
0.1300
0.1866
0.1833
r
5.0
0.1663
0.1689
0.3617
0.3202
0.5004
0.5742
c
5.0
0.0705
0.0704
0.1583
0.1569
0.2250
0.2210Synthetic data assimilation requires the “true” statext(ti), the background (first guess) state xb(ti), and observations y0(ti) inside the assimilation window t0,tN as well as error covariance matrices of the prior guess B0 and observations Ri. The length of data assimilation window should be defined as well. The “true” and background trajectories on the data assimilation interval t0,tN represent some portions of xt(ti) and xb(ti), respectively. Observations should be provided every 5–10 time steps inside the assimilation window and can be generated by adding the Gaussian random noise (with zero mean and variance σ02=0.2) to the “true” state. Since observation grid and model grid are the same, observation operator H is simply an identity mapping. To take into consideration the background covariances, for simplicity, the assumption B0=σb2I can be used, where σb2=0.2 is the variance of background error and I is the identity matrix. Under assumption that the observation quality is the same for all variables, the observation covariance matrix can be defined as Ri=R=σ02I.Testing the TL model and its adjoint is required to ensure the convergence of the minimization algorithm in data assimilation procedures. Ifζδx is a small perturbation of the model state, then(30)Mx+ζδx-Mx≈Mxζδx.To verify the applicability of the TL model on the time interval t0,tN, the relative error (31)eR=Mx+ζδx-MxMxζδxshould be calculated. The TL model is valid if eR→0 when ζ→0. The results of numerical experiments showed that the TL model passed this test with eR tending towards zero (Table 2). The TL adjoint correctness can be tested by verification of the inner product identity(32)Mδx,Mδx=δx,MTMδx.It was found that this equality is essentially correct: the difference was observed only in the 7th digit, which is consistent with a round-off error. The second test to verify the adjoint model is the so-called gradient test [45], which aims to compare a finite difference representation of the gradient of 4D-Var cost function (13) with the gradient obtained via adjoint model ∇Jx0. A linear Taylor approximation of the cost function can be written as (33)Jx0+ζδx≈Jx0+ζδxT∇Jx0.Let us introduce the following function:(34)Ψζ=Jx0+ζδx-Jx0ζδxT∇Jx0.If the gradient is estimated correctly then the function Ψζ→1 as ζ→0. The perturbation vector δx is taken to be [45](35)δx=∇Jx0∇Jx0,where · is the L2 norm. Table 3 manifests the success of the gradient test.Table 2
Results of verification of tangent linear model forc=0.8 and δx=10-2×x0.
ζ
e
R
1
0.9182066027544249
10−1
0.9997279965782743
10−2
0.9999925468155463
10−3
0.9999991929611531
10−4
0.9999999534965012
10−5
0.9999999911217883
10−6
0.9999999914427087
10−7
0.9999997435447022Table 3
Results of verification of 4D-Var cost function gradient forc=0.8 and δx=10-2×x0.
ζ
Ψ
ζ
log
10
⁡
Ψ
ζ
-
1
10−4
0.8727731981461396
−0.8954213897009056
10−5
0.9975483775420343
−2.6105464089686840
10−6
0.9998765512756632
−3.9085133935113787
10−7
0.9999884562844123
−4.9376543818657845
10−8
0.9999979865432855
−5.6960577024708600
10−9
0.9999998912431426
−6.9635433500893210
10−10
0.9999999244103234
−7.1215375125562080


## 6. Sensitivity of the System with respect to Parameters

According to the sensitivity theory [44], general solutions of sensitivity equations for oscillatory nonlinear dynamical systems grow unbounded as time tends to infinity; therefore, sensitivity functions calculated by conventional approaches have a high degree of uncertainty. The reason is that nonlinear dynamical systems that exhibit chaotic behavior are very sensitive to its initial conditions. Thus, the solutions to the linearized Cauchy problem (7) grow exponentially as δxt≈δx0eλt, where λ>0 is the leading Lyapunov exponent. As a result, calculated sensitivity coefficients contain a fairly large error [35–37]. To illustrate this point, let us explore the sensitivity of model output to changes in the coupling strength parameter. Let us introduce the following sensitivity coefficients: (36)S1c=∂x∂c,S2c=∂y∂c,S3c=∂z∂c,S4c=∂X∂c,S5c=∂Y∂c,S6c=∂Z∂c.The corresponding sensitivity equations can be written as(37)S˙1c=σS2c-S1c-cS4c-X,S˙2c=rS1c-S2c-xS3c-zS1c+cS5c+Y,S˙3c=xS2c+yS1c-bS3c+cS6c+Z,S˙4c=εσS5c-S4c-cS1c-x,S˙5c=εrS4c-S5c-XS6c-ZS4c+cS2c+y,S˙6c=ε(XS5c+YS4c-bS6c)-cS3c-z.Sensitivity coefficients can be introduced for any particular model parameter. Since the parameter vector α consists of five components (σ,r,b,c, and ε), five sets of sensitivity equations can be derived from the model equations (2). The dynamics of sensitivity coefficients (36) can be traced by solving the sensitivity equations (37) along with the nonlinear model (2).Sensitivity coefficients (36), calculated on the time interval [0,20], are shown in Figure 4. Envelopes of these coefficients grow over time while sensitivity coefficients themselves exhibit oscillating behavior. Figure 5 provides more detailed information on the evolution of sensitivity coefficients (36) calculated on the time interval [0,5]. It is known that sensitivity coefficient is a measure of the change in state variable due to the variation of the estimated parameter. Unfortunately, obtained sensitivity coefficients are inherently uninformative and misleading. We cannot make a clear conclusion from them about system sensitivity to variations in the parameter c. In this regard, the average values of sensitivity functions ∇αJα over a certain period of time can be considered as one of the most important measures of sensitivity, where J is a generic objective function (20). However, the gradient ∇αJα cannot be correctly estimated within the framework of conventional methods of sensitivity analysis since for chaotic systems it is observed that [35, 36](38)∇αJα≠∇αJα.This is because the integral (39)I=limT→∞∫0Tlimδα→0Jα+δα-Jαδαdtdoes not possess uniform convergence and two limits (T→∞иδα→0) would not commute.Figure 4
Time dynamics of sensitivity functions with respect to parameterc on the time interval [0,20] for c0=0.9.Figure 5
Time dynamics of sensitivity functions with respect to parameterc on the time interval [0,5] for c0=0.9.Similar results were obtained when we considered the influence of variations in the parameterr on the system dynamics. This parameter plays an important role in the formation of system’s dynamical structure and transition to chaotic behavior [32]. Figure 6 shows the differences between components of state vector xr0 obtained with r=r0=28 and xr0+δr obtained with r=r0+δr, where δr=0.01r0=0.28. Even a small perturbation in the parameter r generates a tangible difference between corresponding state variables. Let us define the following sensitivity coefficients:(40)S1r=∂x∂r,S2r=∂y∂r,S3r=∂z∂r,S4r=∂X∂r,S5r=∂Y∂r,S6r=∂Z∂r.The associated system of sensitivity equations can be written as (41)S˙1r=σS2r-S1r-cS4r,S˙2r=x+rS1r-S2r-xS3r+zS1r+cS5r,S˙3r=xS2r+yS1r-bS3r+cS6r,S˙4r=εσS5r-S4r-cS1r,S˙5r=εX+rS4r-S5r-XS6r+ZS4r+cS2r,S˙6r=εXS5r+YS4r-bS6r-cS3r.Figure 6
Difference betweenxr0 and xr0+0.01r0 for fast and slow subsystems.Envelopes of calculated sensitivity coefficients (40) grow over time and sensitivity coefficients demonstrate the oscillating behavior (Figures 7 and 8). Obtained sensitivity coefficients are also uninformative and inconclusive. The “shadowing” approach for estimating the system sensitivity to variations in its parameters [36, 37] allows us to calculate the average sensitivities ∇αJα and therefore to make a clear conclusion with respect to the system sensitivity to its parameters. A detailed description of two variants of the “shadowing” approach is provided in an appendix, and some results of numerical experiments are presented below.Figure 7
Time dynamics of sensitivity functions with respect to parameterr on the interval [0,25] for c0=0.9.Figure 8
Time dynamics of sensitivity functions with respect to parameterr on the time interval [0,5] for c0=0.9.The main problem arising in the “shadowing” method is to calculate the pseudotrajectory. We consider two sets of numerical experiments: weak coupling (c=0.01) and strong coupling (c=0.8) between fast and slow systems. Fast and slow variables that correspond to the original and pseudoorbits are shown in Figures 9 and 10 when the coupling strength parameter c=0.01. The Least Square Shadowing variant of the “shadowing” approach was used to calculate pseudotrajectories. The differences between state variables corresponding to the original and pseudotrajectories of the fast and slow systems are plotted in Figure 11. These figures show that the calculated pseudoorbits are close to corresponding true trajectories over a specified time interval, demonstrating the shadowability. The strong coupling does not introduce significant qualitative and quantitative changes in the behavior of pseudotrajectories with respect to the true orbits. The original and pseudo fast and slow variables for c=0.8 are shown in Figures 12 and 13, and the differences between these state variables are presented in Figure 14. Sensitivity estimates with respect to the parameter r calculated over the time interval [0,20] for different values of coupling strength parameter are shown in Table 4. The most sensitive variables are z and Z. The sensitivity of variables x, y, X, and Y with respect to r is significantly less than variables z and Z.Table 4
Sensitivity estimates of fast and slow variables with respect to parameterr.
c
∂
Z
/
∂
r
∂
Y
/
∂
r
∂
X
/
∂
r
∂
z
/
∂
r
∂
y
/
∂
r
∂
x
/
∂
r
1
.
0
1.10
0.05
0.01
1.08
0.04
0.03
0.8
0.69
0.08
0.03
1.02
0.07
0.07
0.4
0.95
0.03
−0.01
1.03
0.09
0.09
0.15
0.91
−0.08
−0.09
1.01
−0.01
−0.01
1
0
-
4
1.04
−0.02
−0.03
1.02
−0.01
−0.01Figure 9
Original trajectory (in red) and pseudoorbit (in blue) for the fastz and slow Z variables for c0=0.01.Figure 10
Original trajectory (in red) and pseudoorbit (in blue) for fastx and slow X variables for c0=0.01.Figure 11
Difference between variables that correspond to the original trajectory and pseudoorbit forc0=0.01.Figure 12
Original trajectory (in red) and pseudoorbit (in blue) for the fastz and slow Z variables for c0=0.8.Figure 13
Original trajectory (in red) and pseudoorbit (in blue) for fastx and slow X variables for c0=0.8.Figure 14
Difference between variables that correspond to the original trajectory and pseudoorbit forc0=0.8.


## 7. Concluding Remarks

We considered a coupled nonlinear dynamical system, which is composed of fast (the “atmosphere”) and slow (the “ocean”) versions of the well-known Lorenz model. This low-order mathematical tool allows us to mimic the atmosphere-ocean system and therefore serves as a key part of a theoretical and computational framework for the study of various aspects of coupled 4D-Var procedures. Numerical models used to predict the weather are highly nonlinear but tangent linear approximations and their adjoints are used in VDA algorithms. Linear approximation of strongly nonlinear NWP models and also uncertainties in their numerous parameters generate errors in the initial conditions obtained via data assimilation systems. The influence of parameter uncertainties on the results of data assimilation can be studied using sensitivity analysis.We discussed conventional methods of sensitivity analysis and their inefficiency with respect to calculating sensitivity coefficients for chaotic dynamics. To calculate sensitivity coefficients with acceptable accuracy, the sensitivity analysis method [36, 37], developed on the basis of theory of shadowing of pseudoorbits in dynamical systems [38, 39], was applied. Previously, this method was used to analyze the sensitivity of the periodic van der Pol oscillator, the original Lorenz system, and simplified aeroelastic model that exhibit both periodic and chaotic regimes.Calculated sensitivity coefficients obtained via conventional methods and the “shadowing” approach are presented and discussed. It was shown that envelopes of sensitivity coefficients obtained by conventional methods grow over time and the coefficients themselves exhibit the oscillating behaviour. Using the “shadowing” method allows us to calculate the average sensitivity functions (coefficients) that can be easily interpreted.In conclusion, two comments should be highlighted.(1) The shadowing property of dynamical systems is a fundamental feature of hyperbolic systems that was first discovered by Anosov [46] and Bowen [47]. However, most physical systems are nonhyperbolic. Despite the fact that much of shadowing theory has been developed for hyperbolic systems, there is evidence that nonhyperbolic attractors also have the shadowing property (e.g., [48–51]). In theory this property should be verified for each particular dynamical system, but this is more easily said than done.(2) The applicability of the shadowing method for sensitivity analysis of modern atmospheric and climate models is a rather complicated problem since these models are quite complex and they contain numerous input parameters. Thus, further research and computational experiments are required. However, we are confident that, by using the basic ideas of the shadowing method, it is possible to better understand the sensitivity analysis of atmospheric models of various levels of complexity.



---

*Source: 753031-2015-08-06.xml*