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

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.


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].
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][19][20][21][22][23][24][25][26][27][28][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][36][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.

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]): where lower case letters represent the fast subsystem and capital letters the slow subsystem, , , and  are the parameters of L63 model,  is a coupling strength parameter for the  and  variables,   is a coupling strength parameter for ,  is an "uncentering" parameter, and  is a parameter representing the amplitude scale factor.The value  = 1 indicates that two systems have the same amplitude scale.Thus, the state vector of the coupled model (1a) and (1b) is x = (,,,,,)  and the model parameter vector is  = (, , , , ,   , , )  .Without loss of generality, we can assume that  = 1,  = 0, and  =   , then  = (, , , , )  and the system (1a) and (1b) can be rewritten in the operator form: where the nonlinear uncoupled operator L and linear coupled operator Q are represented by the following matrices: The unperturbed parameter values are taken as Chosen values of , , and  correspond to chaotic behaviour of the L63 model.For  = 10 and  = 8/3 the critical value of parameter  is 24.74, which means that any value of  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 parameter  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  tends to zero, then system (2) has two positive, two zero, and two negative Lyapunov exponents.The influence of coupling strength parameter  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 parameter c.At about  ≈ 0.8 they approach the -axis and at about  ≈ 0.95 negative values.Thus, for  > 0.95, the dynamics of both fast and slow subsystems become phase-synchronous [40].For  > 1.0, a limit circle dynamical regime is observed since all six exponents become negative.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 {  } −1 =0 is defined as 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: Signal analysis commonly uses the normalized ACF, defined as () = ()/(0).ACF plots for realizations of dynamic variables  and , and  and  calculated for different values of the coupling strength parameter  are presented in Figures 2 and 3, respectively.For relatively small parameter  ( < 0.4), the ACFs for both  and  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  is less than that of the fast variable .The ACF envelopes for variables  and  also decay almost exponentially from the maximum to zero.For coupling strength parameter on the interval 0.4 <  < 0.6 the ACF of the fast variable  becomes smooth and converges to zero.At the same time, the envelopes of the ACFs of variables , , and  demonstrate a fairly rapid fall, indicating the chaotic behaviour.As the parameter  increases, the ACFs become periodic and their envelopes decay slowly with time, indicating transition to regularity.For  > 0.8 calculated ACFs show periodic signal components.

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: Here x() 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),  is a nonlinear vectorvalued function, such that  : X × P × T → X, and x 0 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  discrete nonlinear dynamical system given by the equation where x  ∈ R  is the -dimensional state vector representing the complete set of the model variables that determine the internal state of the atmospheric model at time   , M ,+1 is the nonlinear operator that indirectly contains model parameters and propagates the state vector from time   to time  +1 for  = 0, . . .,  − 1, and    is the model errors.It is usually assumed that model ( 8) is "perfect" (   = 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 x 0 , (8) uniquely specifies the orbit of the dynamical system.Let y 0  ∈ R  be the -dimensional vector of observations measured at a discrete time   ,  = 0, . . .,  that are linked to the system state via the following equation: where H  : R  → R  is the nonlinear observation operator that maps the state vector to observation space.It is usually assumed that the observation errors  0  are unbiased, serially uncorrelated, and normally distributed with known covariance matrices R  ∈ R × .Suppose that at the initial time the prior (background) model state x  0 is known and represents the "best" estimate of the "true" state x  0 before any observations are taken.This background state is provided by a previous forecast.It is assumed that x  0 has unbiased and normally distributed errors   with known covariance matrix B 0 ∈ R × : Given the observations y   at time   , the corresponding observation error covariance matrices R  ( = 0, . . ., ), the background initial state x  0 , and the error covariance matrix B 0 , the 4D-Var data assimilation seeks to minimize, with respect to x 0 , a certain cost function (x) expressing the "distance" between observations and corresponding model state using the model equations as constraints: subject to x satisfying the set of the "ideal" (   = 0) model equations with initial state x 0 : The 4D-Var cost function is usually written as (e.g., [17,18]): 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: where M  0, is the adjoint of the linearized model operator M 0, = M  0, (x  ) and H   is the adjoint of the linearized observation operator H  = H   (x  ).If the model is not "perfect" then we need to take into account the model errors   , which are sometimes taken as Gaussian noise: where Q is a model error covariance matrix.Thus, we obtain the weakly constrained 4D-Var data assimilation and the following term  Advances in Meteorology 7 should 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.

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 vector   with respect to some model parameter   [30,31,44]: where   is the infinitesimal perturbation of parameter   around some fixed point  0  .Differentiating (8) with respect to , we obtain the set of nonhomogeneous ODEs, the socalled sensitivity equations, which can be written as where S  = (x/  ) = ( 1, ,  2, , . . .,  , )  is the sensitivity vector with respect to parameter   , D  = ( 1 /  ,  2 /  , . . .,   /  )  , and M is a Jacobian matrix.Thus, to analyze the sensitivity of system (7) with respect to parameter   one can solve the following set of differential equations with given initial conditions: 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 time  =  0 by integrating the nonlinear model equations (7) for a long enough time range [ 0 , ], starting from random initial conditions.(2) Solve the nonlinear model equations ( 7) to calculate a trajectory x(),  ∈ [ 0 , ].
(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): 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 x 0 quantifies the influence of parameters on the model output results.In particular, the effect of the th parameter can be estimated as follows: where   is the variation in parameter  0  .Note that

Advances in Meteorology
This approach is acceptable for low-order models.However, the accuracy of sensitivity estimates strongly depends on choice of the perturbation   .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   [30,31].The Gâteaux differential for the objective function (20) has the following form: 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 x 0 (), we obtain the following system of variational equations, the so-called tangent linear model, for calculating x: Then using (24) we can calculate the variation J.Since J(x 0 ,  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  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: where the vector function x * is the solution of adjoint model 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 time  0 by integrating the nonlinear model equations (7) for a long enough time range [ 0 , ], starting from random initial conditions.
(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.

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   = −20 with the initial conditions x (  ) = (0.01, 0.01, 0.01, 0.02, 0.02, 0.02)  (28) and finished at  0 = 0 to guarantee that the calculated model state vector x 0 = 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 x  (  ) and the "true" state x  (  ) at the verification time   are obtained by integrating the model equations over the time interval [ 0 ,   ] with unperturbed parameter vector  0 and initial conditions x 0 .Then, the forecast trajectory x  (  ) and the forecast state x  (  ) at   are obtained by integration of the "forecast" model ( 2) with initial conditions x 0 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: As an example, variations in parameters  and  are considered and forecast is made for two time periods   : 2.5 and 5.0 of nondimensional time units.The forecast error ( 29) is estimated at time   .Table 1 shows the results of forecast verifications.It is obvious that the less the forecast error measure (x  ), 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   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  influences the forecast accuracy almost twice as much as parameter .
Synthetic data assimilation requires the "true" state x  (  ), the background (first guess) state x  (  ), and observations y 0 (  ) inside the assimilation window [ 0 ,   ] as well as error covariance matrices of the prior guess B 0 and observations R  .The length of data assimilation window should be defined as well.The "true" and background trajectories on the data assimilation interval [ 0 ,   ] represent some portions of x  (  ) and x  (  ), 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  2 0 = 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 B 0 =  2  I can be used, where  2  = 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 R  = R =  2 0 I. 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 To verify the applicability of the TL model on the time interval [ 0 ,   ], the relative error should be calculated.The TL model is valid if   → 0 when  → 0. The results of numerical experiments showed that the TL model passed this test with   tending towards zero (Table 2).The TL adjoint correctness can be tested by verification of the inner product identity ⟨Mx, Mx⟩ = ⟨x, M  Mx⟩ .
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 ∇(x 0 ).A linear Taylor approximation of the cost function can be written as Let us introduce the following function: If the gradient is estimated correctly then the function Ψ() → 1 as  → 0. The perturbation vector x is taken to be [45]  where ‖ ⋅ ‖ is the  2 norm.Table 3 manifests the success of the gradient test.

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 ‖x()‖ ≈ ‖x(0)‖  , where  > 0 is the leading Lyapunov exponent.As a result, calculated sensitivity coefficients contain a fairly large error [35][36][37].To illustrate this point, let us explore the sensitivity of model output to changes in the coupling Advances in Meteorology strength parameter.Let us introduce the following sensitivity coefficients: The corresponding sensitivity equations can be written as Sensitivity coefficients can be introduced for any particular model parameter.Since the parameter vector  consists of five components (, , , , 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 .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] This is because the integral does not possess uniform convergence and two limits ( → ∞ K  → 0) would not commute.Similar results were obtained when we considered the influence of variations in the parameter  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 x( 0 ) obtained with  =  0 = 28 and x( 0 + ) obtained with  =  0 + , where  = 0.01 0 = 0.28.Even a small perturbation in the parameter  generates a tangible difference between corresponding state variables.Let us define the following sensitivity coefficients: The associated system of sensitivity equations can be written as 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 ⟨∇  ()⟩ 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.
The main problem arising in the "shadowing" method is to calculate the pseudotrajectory.We consider two sets of numerical experiments: weak coupling ( = 0.01) and strong coupling ( = 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  = 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  = 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  calculated over the time interval [0, 20] for different values of coupling strength parameter are shown in Table 4.The most sensitive variables are  and .The sensitivity of variables , , , and  with respect to  is significantly less than variables  and .

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][49][50][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.

Appendix
The novel sensitivity analysis method for chaotic dynamical systems developed in [36,37] is based on the theory of pseudoorbit shadowing in dynamical systems [38,39], which is one of the most rapidly developing components of the global theory of dynamical systems and classical theory of structural stability [52].Naturally, pseudo-(or approximate-) trajectories arise due to the presence of round-off errors, method errors, and other errors in computer simulation of dynamical systems.Consequently, we will not get an exact trajectory of a system, but we can come very close to an exact solution and the resulting approximate solution will be a pseudotrajectory.The shadowing property (or pseudoorbit tracing property) means that, near an approximate trajectory, there exists the exact trajectory of the system considered, such that it lies uniformly close to a pseudotrajectory.The shadowing theory is well developed for the hyperbolic dynamics, which is characterized by the presence of expanding and contracting directions for derivatives.The study of shadowing problem was originated by Anosov [46] and Bowen [47].
Let (, dist) be a compact metric space and let  :  →  be a homeomorphism (a discrete dynamical system on ).A set of points  = {  :  ∈ Z} is a -pseudotrajectory ( > 0) of  if dist ( +1 ,  (  )) < ,  ∈ Z. (A.1) Here the notation dist(⋅, ⋅) denotes the distance in the phase space between two geometric objects within the brackets.We say that  has the shadowing property if given  > 0 there is  > 0 such that for any -pseudotrajectory  = {  :  ∈ Z} there exists a corresponding trajectory  = {  :  ∈ Z}, which -traces ; that is dist (  ,   ) < ,  ∈ Z. (A. 2) The shadowing lemma for discrete dynamical systems [53] states that, for each  > 0, there exists  > 0 such that each -pseudotrajectory can be -shadowed.
The definition of pseudotrajectory and shadowing lemma for flows (continuous dynamical systems) [38] are more complicated than for discrete dynamical systems.Let Φ  : R ×  →  be a flow of a vector field  on .A function  : R →  is a -pseudotrajectory of the dynamical system Φ  if the inequalities dist (Φ  (,  ()) ,  ( + )) <  (A.3) hold for any  ∈ [−1, 1] and  ∈ R. The "continuous" shadowing lemma ensures that, for the vector field  generating the flow Φ  , the shadowing property holds in a small neighborhood of a compact hyperbolic set for dynamical system Φ  .It is very important to note that the shadowing problem for continuous dynamical systems requires reparameterization of shadowing trajectories.This is the case because for continuous dynamical systems close points of pseudotrajectory and true trajectory do not correspond to the same moments of time.A monotonically increasing homeomorphism ℎ : R → R such that ℎ(0) = 0 is called a reparameterization and denoted by Rep.For  > 0, Rep() is defined as follows [38]: For simplicity, we will consider a generic autonomous dynamical system with one parameter : The new sensitivity analysis method [36,37] is based on the "continuous" shadowing lemma and the following two basic assumptions.
(a) Model state variables are considered over long time interval  ∈ [0, ], where  → ∞, and an averaged With these assumptions, we can use the arbitrarily chosen trajectory of the system to trace the state variables along the orbit and calculate the performance measure ().For example, the arbitrary trajectory () can be chosen as a solution of the model equation, such that it is located nearby a certain reference trajectory   ().Taking into account the shadowing lemma, the closest orbit () to   () satisfies the following constrained minimization problem [37]  where  is the parameter that provides the same order of magnitude of the two terms in the integrand and () is a time transformation.Several numerical algorithms can be used to solve the linearized LSS problem (A.9).One such method is based on variational calculus, which is used to derive optimality conditions representing a system of linear differential equations that are then discretized and solved numerically to calculate variables  and  [37].Let Δ = / be a uniform discretization time step, then denoting  +1/2 =   (( + 1/2)Δ),  −1/2 =   (( − 1/2)Δ),  = 0, . . .,  − 1, and using the trapezoidal rule to approximate the time derivative of  and , we can obtain the following discrete Karush-Kuhn-Tucker (KKT) system [37,54] . . . (A.12) where (A.13) The KKT system can be solved using, for example, iterative methods.Gauss elimination approach for solving the KKT system was considered in [37].The convergence of the LSS method was proved in [55].Calculated variables  and  are then used to determine the sensitivity estimate (A.11).The LSS algorithm is summarized as follows.
(4) Solve the KKT system to obtain variables   and   .
Another method that also uses the concept of pseudoorbit shadowing in dynamical systems for sensitivity analysis of chaotic oscillations was developed in [36].This method In general, the procedure for solving a sensitivity analysis problem represents the following set of steps.
(1) Obtain initial conditions of the system attractor by integrating the model equations (A.5) from   = −20 to  0 = 0, starting from random initial conditions.

Figure 2 :
Figure 2: Autocorrelation functions for dynamic variables  and  for different parameter .

Figure 3 :
Figure 3: Autocorrelation functions for dynamic variables  and  for different parameter .

6 Figure 4 :
Figure 4: Time dynamics of sensitivity functions with respect to parameter  on the time interval [0, 20] for  0 = 0.9.

Figure 5 :
Figure 5: Time dynamics of sensitivity functions with respect to parameter  on the time interval [0, 5] for  0 = 0.9.

Figure 7 :
Figure 7: Time dynamics of sensitivity functions with respect to parameter  on the interval [0, 25] for  0 = 0.9.

Figure 8 :
Figure 8: Time dynamics of sensitivity functions with respect to parameter  on the time interval [0, 5] for  0 = 0.9.

Figure 9 : 20 Figure 10 :
Figure 9: Original trajectory (in red) and pseudoorbit (in blue) for the fast  and slow  variables for  0 = 0.01.

Table 1 :
Relative errors in energy norm (subscript at parameters denotes the nondimensional time unit).

Table 2 :
Results of verification of tangent linear model for  = 0.8 and x = 10 −2 × x 0 .

Table 4 :
Sensitivity estimates of fast and slow variables with respect to parameter .
[37]tisfies the model equation at a different .The desired sensitivity estimate ∇  ⟨() ()⟩ can be computed by solving the following linearized LLS problem[37]:The time-dependent parameter  is called a time-dilation variable and corresponds to the time transformation from the shadowing lemma.Using  and , we can estimate the desired sensitivity (the derivative of the linearized performance measure (A.8) with respect to the parameter ): The second term in the integrand describes reparameterization.Problem (A.7) is called the nonlinear Least Square Shadowing (LSS) problem, and its solution, denoted by () (, ) and  ()  (, ), is a : 1  1  1  2  2  2