Nonlin. Processes Geophys., 24, 922, 2017 www.nonlin-processes-geophys.net/24/9/2017/ doi:10.5194/npg-24-9-2017 Author(s) 2017. CC Attribution 3.0 License.
Zhe An1, Daniel Rey1, Jingxin Ye1, and Henry D. I. Abarbanel1,2
1Department of Physics, University of California, San Diego, 9500 Gilman Drive, La Jolla, CA 92093-0374, USA
2Marine Physical Laboratory (Scripps Institution of Oceanography) University of California, San Diego, 9500 Gilman Drive, La Jolla, CA 92093-0374, USA
Correspondence to: Zhe An ([email protected])
Received: 14 March 2016 Published in Nonlin. Processes Geophys. Discuss.: 22 March 2016 Revised: 1 September 2016 Accepted: 23 September 2016 Published: 16 January 2017
Abstract. The problem of forecasting the behavior of a complex dynamical system through analysis of observational time-series data becomes difcult when the system expresses chaotic behavior and the measurements are sparse, in both space and/or time. Despite the fact that this situation is quite typical across many elds, including numerical weather prediction, the issue of whether the available observations are sufcient for generating successful forecasts is still not well understood. An analysis by Whartenby et al. (2013) found that in the context of the nonlinear shallow water equations on a plane, standard nudging techniques require ob-serving approximately 70 % of the full set of state variables. Here we examine the same system using a method introduced by Rey et al. (2014a), which generalizes standard nudging methods to utilize time delayed measurements. We show that in certain circumstances, it provides a sizable reduction in the number of observations required to construct accurate estimates and high-quality predictions. In particular, we nd that this estimate of 70 % can be reduced to about 33 % using time delays, and even further if Lagrangian drifter locations are also used as measurements.
1 Introduction
The ability to forecast the complex behavior of global circulation in the coupled earth system lies at the core of modern numerical weather prediction (NWP) efforts. Successfully predicting such behavior requires both a good model
Estimating the state of a geophysical system with sparse observations: time delay methods to achieve accurate initial states for prediction
of the underlying physical processes as well as an accurate estimate of the state of the model at the end of the analysis or observation window. When the model is chaotic, even if it is known precisely, the accuracy of the prediction depends on the accuracy of the initial state estimate. This is due to sensitive dependence on the initial conditions, which was rst identied by Lorenz (1963).
Here we consider an idealized situation where a perfect dynamical model describes the deterministic time evolution of a set of D state variables. We assume that L measurements are made at uniform time intervals [Delta1]t within an observation window of length T , so the total number of distinct measurements is L [notdef] (T /[Delta1]t + 1). Our main concern here is the
case where the measurements are sparse in state space, so (L D).
This situation of high-dimensional dynamics and sparse measurements is typical in the process of examining the consistency of observed data and quantitative models of complex nonlinear systems: from functional nervous systems to genetic transcription dynamics, among many other examples (Abarbanel, 2013). Although the methods described here have broad applicability across the quantitative study of the underlying physical or biological properties appearing in many complex systems, our discussion will focus solely on a specic geophysical system: the shallow water equations. As discussed by Cardinali (2013), operational NWP models at the European Centre for Medium-Range Weather Forecasting (ECMWF, 2013) now contain upwards of 108 degrees of freedom. These models are analyzed using 34 [notdef]107 daily
Published by Copernicus Publications on behalf of the European Geosciences Union & the American Geophysical Union.
10 Z. An et al.: Estimating the state of a geophysical system with sparse observations
observations, a large portion of which are often discarded.
Given the scale of these calculations, the question of whether the remaining observations are in fact sufcient for producing accurate analyses and forecasts is of considerable importance.
To clarify these ideas we refer to the observability study given by Whartenby et al. (2013), which evaluated the performance of familiar nudging methods on chaotic, shallow water ow. The ow was simulated on a plane dened by a square grid with uniform spacing N[Delta1] and periodic boundary conditions, and driven by Ekman pumping. Poor predictions were obtained unless the height variable h and at least one of the two velocity variables u, v at each of the N[Delta1] [notdef] N[Delta1]
grid points were measured. In other words, accurate forecasts required direct observation of roughly 70 % of the 3N2[Delta1] dynamical variables.
This lower bound was termed the critical number of measurements Ls required to synchronize the model with the data. Synchronization occurs when the error between the model state estimate and the data drops below a prescribed threshold, which is typically below the magnitude of the observation noise. It depends on a number of factors, including the type of observation network, the signal to noise ratio, properties of the model such as the number and magnitude of its Lyapunov exponents, as well as the choice of algorithm.Strong constraint 4DVar for instance, which is now standard practice in data assimilation (Rabier et al., 2000), encounters serious difculty when the length of the window is long relative to the timescale of the chaos (Pires et al., 1996). In this case, the algorithm will not produce adequate forecasts even with full observations L = D. Despite this, however,
the lowest estimates of Ls appear remarkably consistent between nudging methods and xed interval formulations of 4DVar with both strong and weak constraints (Abarbanel et al., 2009; Abarbanel, 2013; Quinn and Abarbanel, 2010).
Here we examine what can be done when L < Ls. Specically, we will show that, using the method introduced by Rey et al. (2014a, b), which modies a standard nudging technique to include additional information in the time delays of the observations, the estimate of 70 % given by Whartenby et al. (2013) can be reduced to roughly 33 %, and can be even further reduced if positional observations from Lagrangian drifters are also used. These outcomes suggest that time delays may be useful for reducing the number of required observations to meet the practical constraints of operational NWP. However, further testing with more realistic models, observations, and noise is required to verify this claim.
2 Time delayed nudging
We now briey discuss the concept of time delayed nudging; further details can be found in Rey et al. (2014a, b). The system is assumed to be described by a mathematical model whose state is given by a D-dimensional vector x(t). The
model denes a dynamical rule for evolving the x(t) in time, which we assume can be represented as a set of ordinary differential equations (ODEs)
dx(t)dt = F (x(t),t). (1)
If the dynamics of the system are given by partial differential equations (PDEs), such as with uids in an earth systems model, these ODEs may be realized by discretizing the PDEs on a grid. It is worth noting however that a non-trivial amount of discretization error is introduced in this process.
Measurements of the physical system are recorded during an observation window 0 t T = N [Delta1]t, where L ob
servations y(tn) are taken at each time tn = n[Delta1]t for n =
0,1,...,N. The measurements y(t) are related to the state vector x(t) through a measurement operator h(t). For simplicity, it is taken here to be a constant, L [notdef] D projection
matrix h(t) = H, so that y(t) = H [notdef] x(t) + noise. The to
tal number of measurements in the observation window is L [notdef] (N + 1).
The overall objective is to estimate the full model state x(T ) at the end of the assimilation window using information from observations, and then use this estimate to predict the systems subsequent behavior for t > T using Eq. (1). The accuracy of these predictions, when compared with additional measured data in the prediction window t > T , serves as a metric to validate both the model and the assimilation method, through which the unobserved states of the system are determined. This establishes a necessary condition on L that is required to synchronize the model output with the data and thereby obtain accurate estimates for the unobserved states of the system.
When the model is known precisely, a familiar strategy for transferring information from the measurements to the model involves the addition of a coupling or control or nudging term to Eq. (1),
dx(t)dt = F (x(t),t) + H [notdef] G(t) [notdef] y(t) H [notdef] x(t)
. (2)
where H denotes the transpose, and G(t) is an L[notdef]L matrix
that is nonzero only at times tn where measurements occur.
For simplicity, when G(t) is non-zero, it is assumed to be constant and diagonal, so the coupling terms only affect the measured states.
This long-standing procedure, known as nudging in the geophysics and meteorology literature, has been shown to fail when the number of measurements at a given time falls below a critical value Ls (Abarbanel et al., 2009). This can be understood by noting that the coupling term perturbs the observed model states, driving them towards the data. With enough observations L, and a sufciently strong coupling G(t), this control term alters the Jacobian of the dynamical system Eq. (2) so that all its (conditional) Lyapunov exponents are negative (Pecora and Carroll, 1990; Abarbanel,
Nonlin. Processes Geophys., 24, 922, 2017 www.nonlin-processes-geophys.net/24/9/2017/
Z. An et al.: Estimating the state of a geophysical system with sparse observations 11
1996; Kantz and Schreiber, 2004). That is, the log of the maximum eigenvalue of the matrix [[Phi1](T ,t0) [notdef][Phi1](T ,t0)]1/2T
is negative, where [Phi1](t,t[prime]) = @x(t)/@x(t[prime]) is the linearized
state transition matrix or tangent linear model. Its time evolution is described by the variational equation
d[Phi1](t,t[prime])
dt = D
eF (x(t),t) [notdef] [Phi1](t,t[prime]) [Phi1]a b(t,t[prime]) = a b, (3)
along the trajectory given by Eq. (2) and
D
eF = DF (x(t),t) H [notdef] G(t) [notdef] His its Jacobian. That is, DFa b(x(t),t) = @Fa(x(t),t)/@xb(t)
and a b is the Kronecker delta, so [Phi1](t[prime],t[prime]) is an identity matrix. This establishes a necessary condition on L required to synchronize the model with the data. Numerical experiments have shown that when this condition is not met, estimates are not accurate and predictions are unreliable (Abarbanel et al., 2009; Abarbanel, 2013). An example of this will be given later in our discussion of geophysical shallow water ow.
It is therefore important to understand, for a given problem, whether L > Ls. If this condition is not satised and additional measurements cannot be made, then we must nd another means to overcome this decit in L. One way to proceed involves the recognition that additional information resides in the temporal derivatives of the observations. In practice, however, this derivative information cannot be measured directly, although it can be approximated via nite differences, for instance by approximating dy(tn)/dt with (y(tn + ) y(tn))/ where is some multiple of the time
differences between measurements. The drawback here is that the derivative operation acts as a high-pass lter, and is thus quite susceptible to noise in the measurements. Alternatively, it has been known for some time in the nonlinear dynamics literature that this additional information in the derivative is also available in the time delay of the measurements, y(tn +). This process can be repeated as many times
as needed to form a DM-dimensional vector of time delays, which we call S(t).
This idea provides the basis for the well-established technique in the analysis of nonlinear dynamical systems, where this structure is employed as a means of reconstructing unambiguous orbits of a partially observable system (see, e.g., Aeyels, 1981a, b; Ma, 1981; Sauer et al., 1991; Takens, 1981; Kantz and Schreiber, 2004; Abarbanel, 1996). By mapping to a proxy space of time delayed observations, one is able to invert the projection associated with measuring L < D components of the underlying dynamics, by using the fact that new information beyond y(tn) lies in y(tn + ). The
derivative operation is just another (albeit less numerically robust) way of accessing this information.
Note that the time delay and the embedding dimension DM are parameters that need to be chosen appropriately for the system, although a number of useful heuristics are available (Abarbanel, 1996). Moreover, Takens (1981) proved that
taking DM > 2DA, where DA is the fractal dimension of the attractor, is sufcient to unambiguously reconstruct the topology of the attractor. It is worth noting however that this condition is only sufcient, and the procedure often succeeds with a considerably smaller value of DM.
In the estimation context, the time delays are used in a slightly different way. Instead of reconstructing the topology of the attractor, they are used to control local instabilities in the dynamics, which cause errors in the analysis to grow. In other words, DM does not need to embed the entire space. Rather, it only needs to be large enough to effectively increase the amount of information transferred from the L measurements to a value above the critical threshold, Ls.
Using this idea Rey et al. (2014a, b) proposed a technique to extract additional information from time delayed observations by constructing an extended state space S(t), created from an L [notdef] DM dimensional vector of the measurements and
its time delays. The components of this time delayed observation vector are denoted by
(4)
where DM is the dimension of the time delayed vector Y (tn), and is the delay, which here is assumed to be a positive integer multiple of [Delta1]t. Also, note that the term delay here is not used in its usual sense. Rather, this method involves a time advanced formulation, which for positive incorporates observations at later times. Both formulations are acceptable however.
The corresponding time delay model vectors S(x(t)) are given by
S(x(t)) =[braceleftBig]
[H [notdef] x(t)],[H [notdef] x(t + )],...,
[H [notdef] x(t + (Dm 1))][bracerightBig]
, (5)
where the values x(t[prime] > t) are constructed by integrating the uncoupled dynamics, Eq. (1), forward in time. The time evolution for S(x(t)) is given by the chain rule,
dS(x(t))dt = DS(x(t)) [notdef] F (x(t),t), (6)
where the Jacobian DS(x(t)) = @S(x(t))/@x(t) with re
spect to x(t) can be computed using the variational Eq. (3), by substituting the Jacobian of the uncoupled model D
eF !
DF . Furthermore, in analogy with Eq. (2), we introduce a control term g(t) in time delay space:
dS(x(t))dt = DS(x(t)) [notdef] F (x(t),t) + g(t) [notdef] Y (t) S(x(t))
. (7)
We then transform back to physical space, by multiplying both sides of this equation by [DS(x(t))]1, to get
dx(t)dt = F (x(t),t) + G(t) [notdef]
DS(x(t))
1
. (8)
www.nonlin-processes-geophys.net/24/9/2017/ Nonlin. Processes Geophys., 24, 922, 2017
[notdef] g(t) [notdef] Y (t) S(x(t))
12 Z. An et al.: Estimating the state of a geophysical system with sparse observations
Note that there are now two control terms, G(t) and g(t), which act in physical and time delay space, respectively.Also, since DS(x(t)) is a (L[notdef]Dm)[notdef]D matrix, it is generally
not square, so its pseudoinverse [DS(x(t))]+ is used.
At each step of the integration of the controlled (nudged) dynamical Eq. (8), the control term perturbs the full state vector in time delay space S(x(t)) toward the time delay measurement vector Y (t), allowing it to extract additional information from the waveform of the existing measurements. The value of this statement will become clearer later on.
Furthermore, in the limit DM = 1 the time delay formula
tion Eq. (8) reduces to the standard nudging control Eq. (2).
Two important differences however are realized when DM >
1. First, information from the time delays of the observations is presented to the physical model equations. And second, all components of the model state x(t) are inuenced by the control term, not just the observed components. This, for example, allows xed parameters p of the model to be estimated as a natural result of the synchronization process by including them as additional state variables, satisfying dp(t)/dt = 0.
Time delay nudging shares considerable overlap with incremental formulations of strong constraint 4DVar (Lewis and Derber, 1985; Talagrand and Courtier, 1987; Courtier et al., 1994). For instance, both methods use a sliding window of observations and compute the control (nudging) perturbation by minimizing the squared magnitude of time-distributed innovations [notdef]Y (t)S(x(t))[notdef]2. The use of time ad
vanced observations is also standard practice in strong constraint 4DVar, and is motivated by the fact that the necessary conditions for synchronization require one to control the propagation of errors on the unstable manifold (Trevisan et al., 2010; Palatella et al., 2013). Since these errors are locally described by Eq. (3) as the system evolves forward in time, the time advanced construction is a natural choice.
The main differences between the two methods are as follows.
1. Strong constraint 4DVar does not include the notion of a time delay or embedding dimension.
2. With the time delay method, the analysis is propagated in small increments dt between analyses, and observations are re-used.
3. Time delay nudging uses truncated singular-value decomposition to regularize the solution, while strong constraint 4DVar uses a background term to perform Tikhonov regularization.
Regarding the second point, near the end of the observation window one must either switch to a time delayed formulation or reduce DM appropriately. Here however for simplicity, the end of the observation window is taken so that the last observation y(T + (DM 1)) is always available. Also, the third
point prevents time delay nudging from being applied to sys-
tems of the size used in operational NWP. However, a simplied variation of time delay nudging was recently given by Pazo et al. (2016). This method avoids the variational Eq. (3) and the generalized inverse altogether, and thus requires considerably less computational effort than either time delay nudging or strong constraint 4DVar. It is worth investigating whether this method is also capable of achieving the same reduction in Ls, shown here for geophysical ows.
Furthermore, while we are currently working on unifying the motivating ideas behind time delay nudging with the variational action principle of weak constraint 4DVar, these and other related connections to 4DVar will be given in a subsequent paper. For the moment however, no additional theory will be introduced. Instead, we focus on its application to a core geophysical model: the shallow water equations.
3 Twin experiments
We test our time delay nudging procedure through a series of twin experiments (Durand et al., 2002; Blum et al., 2009; Blum, 2010). After solving the original dynamical Eq. (1) forward from preselected initial conditions x(0), the observation process is simulated by applying the observation operator H to project the state down to the L observed components. Gaussian noise N(0,) is then added to each component to simulate observation error. The estimation process continues as described above until the time t = T . At this
point, the coupling terms g(t) and G(t) are set to zero and the uncoupled dynamics Eq. (1) are integrated forward from the estimated x(T ) to construct a forecast for t > T . Comparing this forecast against additional observations y(t > T ) then tests whether the unobserved states are also accurately estimated.
To simulate the conditions of a true experiment we monitor our progress by calculating the observable synchronization error, namely the root mean square deviation between the data and the observed model states:
SE(tn) =
[radicalbigg]
1 L
[vextendsingle][vextendsingle]H
[notdef] xs(tn) ys(tn) [vextendsingle][vextendsingle]2.
(9)
In this expression, scaled variables have been introduced such that xs[lscript](t) =
x[lscript](t) xmin[lscript](t)
[bracketrightbig]
/
xmax[lscript](t) xmin[lscript](t)
[bracketrightbig]
and
xmin/max[lscript](t) are the minimum or maximum values of x[lscript](t) over the entire assimilation window. The same denition holds for ys[lscript](t). This rescales all data and observed model states to lie in the interval [0,1], so that each state compo
nents contribution to the synchronization error is roughly equal. While this could make the result sensitive to outliers in the data, it did not appear to be an issue here.
It was previously shown by Whartenby et al. (2013) that when the synchronization error Eq. (9) decreases to very small values, the full state x(T ) is accurately estimated and the forecast is quite good. Conversely, when this fails to occur, the full state x(T ) is not well estimated and the predic-
Nonlin. Processes Geophys., 24, 922, 2017 www.nonlin-processes-geophys.net/24/9/2017/
Z. An et al.: Estimating the state of a geophysical system with sparse observations 13
tion is unreliable. In Rey et al. (2014a, b), this contraction of the synchronization error was only observed when the number of time delayed observations L[notdef]DM, and the magnitude
of the coupling matrices g(t), G(t) were large enough.
The precise meaning of this statement will become apparent shortly.
4 Nonlinear shallow water equations
We now describe the application of time delay nudging to a nonlinear model of shallow water ow on a mid-latitude plane. This geophysical uid dynamical model was previously examined by Pedlosky (1987) and Whartenby et al. (2013) as well as many others, and is at the core of earth system ows used in NWP. Of course, operational models contain considerably more detail than this example, and those models also describe the dynamics over a sphere. While we suspect the results presented here for this simplied model will be applicable to more realistic models as well, additional experiments are needed to validate this claim.
As the depth of the coupled atmosphere ocean uid layer (1015 km) is markedly less than the earths radius (6400 km), the shallow water equations for two-dimensional ow provide a good approximation to the uid dynamics of the ocean. Three elds on a mid-latitude plane describe the uid ow [notdef]u(r,t),v(r,t),h(r,t)[notdef]: the northsouth velocity
v(r,t), the eastwest velocity u(r,t), and the height of the uid h(r,t), with r = [notdef]x,y[notdef]. The uid is taken as a single,
constant density layer and is driven by wind stress (r,t) at the surface z = h(r,t) through an Ekman layer. These physi
cal processes satisfy the following dynamical equations with u(r,t) = [notdef]u(r,t),v(r,t)[notdef],@u(r,t)
@t = u(r,t) [notdef] ru(r,t) g rh(r,t) + u(r,t)
(f (y) z) + Ar2u(r,t) [epsilon1] u(r,t), @h(r,t)
@t = r [notdef]
h(r,t)u(r,t)
[bracketrightbig]
z [notdef] curl
2
[bracketleftbig]
cos(! (r(i,j)) + ) + cos(! (r(i,j)) + ) [bracketrightbig]
+H0,
(r,t) f (y)
. (10)
The Coriolis force is linearized about the Equator f (y) =
f0 + y and the wind-stress prole is selected to be (r,t) = [notdef][F/] cos(2 y),0[notdef]. The parameter A represents the viscos
ity in the shallow water layer, [epsilon1] is Rayleigh friction and z is
the unit vector in the z direction. The values we have used for the model parameters are given in Table 1. With these xed parameters the shallow water ow is chaotic, and the largest Lyapunov exponent for this ow is estimated to be max = 0.0325/h 1/31 h by measuring the average growth
rate of random perturbations. The details of this calculation are given by Whartenby et al. (2013).
We have analyzed this ow using the enstrophy conserving discretization scheme given by Sadourny (1975) on a grid of size N2[Delta1] for increasing resolution N[Delta1] = [notdef]16,32,64[notdef].
The total domain size is constant 800 [notdef] 800 km and pe-
riodic boundary conditions are enforced. Using the twin-experiment framework, with simple nudging given in Eq. (2) and a static, uniform observation operator, approximately 70 % of the D = 3N2[Delta1] degrees of freedom were required to
be observed to synchronize the model output with the data (Whartenby et al., 2013). In other words, the height eld and at least one of the velocity elds at each grid point needed to be observed.
Since these results were roughly consistent among the three resolutions tested, we restrict our discussion here to the case where N[Delta1] = 16. Representative plots of the height and
velocity elds are shown in Fig. 1. The total number of degrees of freedom is D = 3N2[Delta1] = 768, for which Whartenby
et al. (2013) estimated Ls 524 = 0.68D. Based on the dis
cussion above and the lectures by Cardinali (2013), we see that this requirement, which is expected to be even higher in practice, exceeds the number of measurements now available by at least a factor of 2.
5 Results with time delay nudging for the shallow water equations
We now demonstrate that the time delay method is capable of reducing Ls, by showing that it can construct successful estimates and predictions without directly observing the horizontal velocity elds. This strategy was shown to fail by Whartenby et al. (2013) with static (DM = 1) nudging.
Thus, we assume height measurements alone are made at each grid point (i,j) for i,j = [notdef]1,2,...,16 = N[Delta1][notdef], so L =
256 < 524 Ls, as estimated by Whartenby et al. (2013).
The initial state x(t0) for the model and the data are taken to have the forms
h(i,j)(t0) = [parenleftbigg]
A0 N[Delta1] [Delta1]Y
u(i,j)(t0) = A0
@ (r(i,j))
@y ,
v(i,j)(t0) = A0
@ (r(i,j))
@x , (11)
with the parameters H0 = 5100 m, A0 = 106 and
(r) = cos
![prime] (r) + [prime][parenrightBig]sin ![prime] (r) + [prime]
[parenrightbig]
. (12)
The functions and , respectively, evaluate the latitude and longitude at the point r(i,j) on the grid. All elds as well as
the variables [notdef]x,y,t[notdef] were scaled by the values in Table (1),
to make them dimensionless. The parameters !, ! , ![prime], ![prime] and , , [prime], [prime] are chosen arbitrarily, so that the phase and period of the initial condition are different for truth and the estimate. Also, although the method is capable of estimating the static model parameters, here they are considered known.
www.nonlin-processes-geophys.net/24/9/2017/ Nonlin. Processes Geophys., 24, 922, 2017
14 Z. An et al.: Estimating the state of a geophysical system with sparse observations
Table 1. Parameters used in the generation of the shallow water data for the twin experiment. All elds as well as [notdef]x,y,t[notdef] were scaled by
the values in the table, so all calculations were done with dimensionless variables.
Parameter Physical quantity Value in twin experiments
[Delta1]t Time step 36 s [Delta1]X Eastwest grid spacing 50 km [Delta1]Y Northsouth grid spacing 50 km H0 Equilibrium depth 5.1 km '0 Central latitude of the plane 3.6
f0 Central value of the Coriolis parameter 5 [notdef] 10
Meridional derivative of the Coriolis parameter 2.0 [notdef] 10
F/ Wind stress 0.2 m2 s3A Effective viscosity 104 m2 s1 [epsilon1] Rayleigh friction 2 [notdef] 10
8 s1
Figure 1. Snapshots of shallow water ow on a 16 [notdef] 16 grid spanning 400 km on each side. Heights and velocity elds are shown: left panel
at the initial time, center panel after 30 min, and right panel after 30 h.
The coupling matrix G(t) is taken to be diagonal, with different weights for the heights and for the velocities. In particular, Gu,v [Delta1]t = 0.5 and Gh [Delta1]t = 1.5 with [Delta1]t = 0.01h =
36 s. The values of Gh are larger than Gu, Gv, since the average height is 5000[notdef]30 m, 3 orders of magnitude higher than
the average velocity 0 [notdef] 5 ms1. The time delay space cou
pling g(t) is taken to be the identity matrix, as all the height measurements are assumed to be known with equal temporal precision throughout the observation window.
The time delay was selected to be = 10[Delta1]t = 0.1 h, in or
der to maintain a balance between numerical stability and the common criterion of independence between the components of S(x(t)). The rst minimum of the average mutual information was also calculated to be 30[Delta1]t using the method
given by Abarbanel (1996). This is reasonably close to our choice, and the results did not change if its value was shifted by a few [Delta1]t.
5.1 Choosing DM
The state was estimated by integrating the coupled differential Eq. (8) from t = 0 to T = 5h = 500[Delta1]t with various
DM = [notdef]1,6,8,10[notdef]. The coupling terms were then switched
off at t = T to generate predictions until t = 500 h.
5 s1
11 m1 s1
Short- and long-term synchronization error (Eq. 9) trajectories SE(t) are plotted in Fig. 2 for various DM. Choosing DM = [notdef]1,6[notdef] yields a synchronization error that remains
around its initial value of 0.005 until the end of the 5 h observation window. After the coupling is switched off, the error rises very rapidly until stabilizing around 0.1 for the remainder of the prediction window. By contrast, for DM =
{8,10[notdef] the synchronization error falls steeply to order 106
within the observation window. It then subsequently rises as exp[ max(t T )], where max 1/31 h agrees with the
largest Lyapunov exponent calculated for this ow. This exponential rate of growth is particularly evident in the long trajectory displayed in the right panel of Fig. 2.
Since DM 8 produces error values several orders of
magnitude smaller than those obtained with DM 6, we ex
pect the state estimates x(T ) obtained with DM 8 to be
quite accurate when compared with the estimates for DM
6. These estimates are now evaluated as they would be in a true experiment, by comparing predictions on the observed heights with additional data. In Fig. 3 the known (black), estimated (red), and predicted (blue) height trajectories are shown for an arbitrarily selected grid point h(6,4)(t). Short-
and long-term prediction trajectories computed with DM = 6
are displayed in Fig. 3s upper panels, respectively. Corre-
Nonlin. Processes Geophys., 24, 922, 2017 www.nonlin-processes-geophys.net/24/9/2017/
Z. An et al.: Estimating the state of a geophysical system with sparse observations 15
Figure 2. Synchronization error SE(t), dened in Eq. (9), computed with DM = [notdef]1,6,8,10[notdef], Gh[Delta1]t = 1.5, Gu[Delta1]t = gv[Delta1]t = 0.5, and =
10[Delta1]t = 0.1 h. Assimilation is performed for t 5 h. Left panel: the couplings are then switched off and predictions are generated using the
original dynamical Eq. (10) until t = 100 h. In the prediction window (t 5), the error in the trajectories grows roughly with the largest
Lyapunov exponent of the system max 1/31 h. Synchronization is evident when DM = [notdef]8,10[notdef] and not for DM = [notdef]1,6[notdef], suggesting that
accurate predictions will be obtained for DM = [notdef]8,10[notdef]. Right panel: the same calculation, but extended to t = 500 h.
Figure 3. Upper left panel: known (black), estimated (red) and predicted (blue) values for the observed height values h(6,4)(t) at the grid point (6,4) for DM = 6. Observations are for 0 t 5 h. Predictions are for 5 t 100 h. Upper right panel: the same calculation for
DM = 6 for a longer prediction window 5 t 500 h. Lower left panel: the same calculation except DM = 8. The prediction window is
5 t 100 h. Lower right panel: the same calculation except DM = 8. The prediction window is 5 t 500 h.
sponding results for DM = 8 are shown in the lower panels.
As anticipated, the predictions for DM = 8 are clearly su
perior to those obtained with DM = 6. In addition, with the
choice DM = 7, some initial conditions synchronized, while
others did not. Further analysis of this boundary case is an interesting area for future study.
The failure of predictions obtained with DM = 6 is a re
sult of poor estimates of the unobserved states (i.e., uid velocities) at t = T . Although in an actual experiment we
would not be able to verify this statement directly, we may do so here. Velocity proles u(6,4)(t) displaying short and long time comparisons between the known (black), estimated
www.nonlin-processes-geophys.net/24/9/2017/ Nonlin. Processes Geophys., 24, 922, 2017
16 Z. An et al.: Estimating the state of a geophysical system with sparse observations
Figure 4. Upper left panel: known (black), estimated (red) and predicted (blue) for the observed x-velocity values u(6,4)(t) at grid point (6,4) for DM = 6. Observations are for 0 t 5 h. Predictions are for 5 t 100 h. Upper right panel: the same calculation for DM = 6 for a
longer prediction window 5 t 500 h. Lower left panel: the same calculation except DM = 8. The prediction window is 5 t 100 h.
Lower right panel: the same calculation except DM = 8. The prediction window is 5 t 500 h.
(red), and predicted (blue) values are given in Fig. 4 for DM = 6 in the upper panels, and for DM = 8 in the lower
panels. We nd the situation is indeed as anticipated; the estimates and predictions are quite unacceptable for DM = 6,
whereas for DM = 8 they are highly accurate. The same im
provement in predictive accuracy was obtained for the other velocity component v(6,4)(t). These results are plotted in Fig. 5.
Predictions were also calculated for DM = 1 and DM =
10, but these results are not shown. They agree with the synchronization error calculations in Fig. 2, in that the predictions generated with DM = 10 are just as accurate as those
for DM = 8. Likewise, predictions with DM = 1 (i.e., simple
nudging) are very poor, in accordance with Whartenby et al. (2013).
5.2 Reducing the coupling strength
In the previous discussion it was suggested that reducing the coupling strength would have a detrimental effect on the quality of the estimation procedure and the resulting prediction. We investigate this now by performing the same calculations as above with DM = 10 but reducing the coupling on
the height so that we have Gh [Delta1]t = Gu [Delta1]t = Gv [Delta1]t = 0.5.
The synchronization error SE(t), shown in Fig. 6, upper left panel, stabilizes to a level 3 orders of magnitude higher than was achieved with Gh [Delta1]t = 1.5, suggesting failure. This is
conrmed in the remainder of Fig. 5, which displays the known (black), estimated (red), and predicted (blue) values for h(6,4)(t), u(6,4)(t), and v(6,4)(t), respectively. Although the height estimate is rather good and the prediction is not terrible, at least for the rst 1520 h after the end of the assimilation window, the unobserved states are clearly not well estimated at any time t T .
This result demonstrates that proper choice of coupling is required. However, we have not developed a systematic way of choosing these values, and it is known from classical results on synchronization that the optimal choice depends on the number and distribution of observations. Furthermore, the fact that the height estimates appear to be rather accurate also emphasizes the point that, in a true experiment, the success of the assimilation procedure must be evaluated against the forecasts not the analyses.
Nonlin. Processes Geophys., 24, 922, 2017 www.nonlin-processes-geophys.net/24/9/2017/
Z. An et al.: Estimating the state of a geophysical system with sparse observations 17
Figure 5. Upper left panel: known (black), estimated (red), and predicted (blue) for the observed y-velocity values v(6,4)(t) at grid point (6,4) for DM = 6. Observations are for 0 t 5 h. Predictions are for 5 t 100 h. Upper right panel: the same calculation for DM = 6 for
a longer prediction window 5 t 500 h. Lower left panel: the same calculation except DM = 8. The prediction window is 5 t 100 h.
Lower right panel: the same calculation except DM = 8. The prediction window is 5 t 500 h.
5.3 Further reducing the number of measurements
In addition, until now we have conveniently chosen to ob-serve the height eld at all L = N2 = 256 grid locations. We
now attempt to reduce L even further, by repeating the analysis with L = 252 and L = 248 height measurements, chosen
at arbitrary grid points. From the results displayed in the upper left panel of Fig. 7, it is evident that for L = 252 rapid and
accurate synchronization is still achieved, while for L = 248
it is not. In addition, the known (black), estimated (red), and predicted values (blue) for h(6,4)(t) are shown in the other panels of Fig. 7 for L = 248 and L = 252, respectively. Re
sults for the unobserved velocity elds agree as well, though
these results are not shown.
Thus, even with time delays, it may not be possible to signicantly reduce the number of required height measurements. We remark however that the overall space of parameters appearing in our study has not been thoroughly explored. Additional renement of the parameters G(t), g(t), DM, and may further reduce this constraint, for instance, by allowing G(t) to be non-diagonal.
5.4 Noise in the observations
We now repeat the above calculations for L = 252 with
Gaussian noise N(0,) added to the height observations.
A comparison is shown in Fig. 8 for = [notdef]0.2,0.5[notdef] and
DM = [notdef]8,10[notdef]. The synchronization error still falls rapidly
within the observation window, although not to O(105), as
in the noiseless case. In the prediction window, it rises in an exponential manner as expected. Furthermore, results fail to synchronize when the magnitude of the noise gets large enough. This transition occurs at roughly = [notdef]1.3,2.0[notdef] for
DM = [notdef]8,10[notdef], respectively. These results were included to
show that the method appears to be relatively robust to small errors in the observations. A more thorough examination of the impact of imperfect observations will be given elsewhere.
5.5 Using drifter data
Another quite important source of observations about ocean ows is being provided by position measurements r(t) of Lagrangian drifters (Mariano et al., 2002). Such observations have been shown to be a good supplement to the traditional observations made on a xed grid (Kuznetsov et al.,
www.nonlin-processes-geophys.net/24/9/2017/ Nonlin. Processes Geophys., 24, 922, 2017
18 Z. An et al.: Estimating the state of a geophysical system with sparse observations
(6,4)(t) at location (6,4), gh [Delta1]t = gu [Delta1]t =
gv [Delta1]t = 0.5. All other parameters are the same. Upper left panel: SE(t) for 0 t 200 h. Upper right panel: known (black), estimated (red),
and predicted (blue) for the observed height values h(6,4)(t) at grid point (6,4) for DM = 10. Observations are for 0 t 5 h. Predictions are
for 5 t 100 h. Lower left panel: known (black), estimated (red), and predicted (blue) for the observed x-velocity values u
Figure 6. Data assimilation results with DM = 10 and reduced coupling on the height component h
(6,4)(t) at grid
point (6,4) for DM = 6. Observations are for 0 t 5 h. Predictions are for 5 t 100 h. Lower right panel: known (black), estimated
(red), and predicted (blue) for the observed y-velocity values v(6,4)(t) at grid point (6,4) for DM = 6. Observations are for 0 t 5 h.
Predictions are for 5 t 100 h.
2003) and they can also be used to estimate an Eulerian velocity eld (Molcard et al., 2003; Piterbarg, 2008; Salman et al., 2006). In this section, we combine the time delay method with a data set from drifter measurements to show that they can provide accurate estimates for the grid state variables
h(r(i,j),t),u(r(i,j),t),v(r(i,j),t)
, without much
additional effort.
We monitor the positions of ND drifters deployed at randomly chosen grid locations and afterwards allowed to move freely to provide spatially continuous measurements between grid points. Drifters were also deactivated when they reached the boundary of the grid, so the number of operational drifters decreases throughout the estimation window. The dynamics of drifters are approximated as two-dimensional uid parcel motion near the surface of the water layer, which are determined by the Lagrangian equations
dr(n)(t)
dt = u(r(n)(t),t), (13)
where r(n)(t) is the position of the nth drifter and this equation was simulated by linear interpolation of the discrete velocity elds (Press et al., 2007; Thomson and Emery, 2014). Hybrid measurements are incorporated into the time delay nudging method by combining the grid variables and the collective drifter positions
(t) = [braceleftBig][bracketleftBig]
r(1)(t)
i,
hr(2)(t)
i,...,
hr(ND)(t)
i [bracerightBig]
into a single hybrid state vector. The corresponding time delayed vectors are Y drifter(t) = [notdef]Y grid(t),Y drifter(t)[notdef] and
Sdrifter(t) = [notdef]Sgrid(t),Sdrifter(t)[notdef], respectively, where
Y drifter(t)
= ndata(t),data(t + ),...,data (t + (DM 1))[bracerightBig] ,
Sdrifter(t)
= nmodel(t),model(t + ),...,model (t + (DM 1))[bracerightBig] .
Nonlin. Processes Geophys., 24, 922, 2017 www.nonlin-processes-geophys.net/24/9/2017/
Z. An et al.: Estimating the state of a geophysical system with sparse observations 19
the synchronization error for L = 208 height observations
and ND = 20 drifter observations, and we show (in blue)
the same synchronization error when L = 208 and ND = 0
drifters are deployed. With L = 208, namely, observing 27 %
of the heights and 20 drifters, the synchronization error converges to a small value within the 5 h observation window. Without drifters, the estimation fails. Furthermore, by increasing the number of drifters to ND = 64 within a 30 min
observation window, synchronization can be achieved with
L = 128 height observations. Snapshots of the elds at dif
ferent times throughout the estimation and prediction window are shown in Fig. 11 for comparison.
Although we have not yet explored how to balance the number of drifters and the number of height (or other) measurements, these preliminary results suggest that positional data from drifters can be useful for improving the observ-ability of the system. In contrast to other approaches in which the drifter data are used to directly interpolate the grid variables (Kuznetsov et al., 2003), our method transfers positional information from the drifters to the estimate through the dynamical model. Whether this approach is valid for real
www.nonlin-processes-geophys.net/24/9/2017/ Nonlin. Processes Geophys., 24, 922, 2017
Figure 7. Synchronization error and known, estimated, and predicted height values for L = 248 height measurements at each observation
time and for L = 252 height measurements at each observation time. Upper left panel: SE(t) for L = 248 and L = 252 over 0 t 5 h in
the observation window, and 5 t 500 h after the couplings are removed. Upper right panel: known (black), estimated (red), and predicted
(blue) values of the height h(6,4)(t) at grid point (6,4) for 0 t 100 h for L = 248. Lower panel: known (black), estimated (red), and
predicted (blue) values of the height h(6,4)(t) at grid point (6,4) for 0 t 100 h for L = 252. This shows the rather sharp transition
between bad predictions (L = 248) and good predictions (L = 252).
We consider cases with ND = 0 (no drifters), as well as
ND = 20 and ND = 64, in addition to L = 208 and L = 128
grid observations. Example plots showing the initial locations for ND = 20 and ND = 64 are given in Fig. 9. While
many runs were successful using the same setup described above, the results were somewhat dependent on where the drifters were initialized. These results are not shown.
The consistency of the results improved by choosing the initial estimate to only magnitude from the true solution, rather than in both phase and frequency as was done above. Specically, for the results reported below the initial conditions of the data data(r(i,j),t0) and hdata(r(i,j),t0) and of the model model(r(i,j),t0) and hmodel(r(i,j),t0) are related by
data(r(i,j),t0) = C0 model(r(i,j),t0) and hdata(r(i,j),t0) =
C0 hmodel(r(i,j),t0). We choose C0 = 1.0 + 0.1 , with se
lected from a uniform distribution in the interval [1,1]. The
velocity elds are found as above, using (r,t0) as a stream function.
In Fig. 10, we show the synchronization error of observed quantities when DM = 8, keeping all other parameters the
same as in the previous calculations. We present (in red)
20 Z. An et al.: Estimating the state of a geophysical system with sparse observations
Figure 8. The effect of noise levels in the initial condition for the solution of the model Eq. (10) on SE(t). We show the results for DM = 8 and 10 for added Gaussian noise N(0,) with = 0.2 and
0.5. For this range of noise levels added to the initial condition for generating the data in our twin experiments, we see that the detailed values of SE(t) change. In the case of both DM = 8 and DM =
10, SE(t) still becomes quite small in the observation window 0
t 5 h, suggesting that predictions for t 5 will remain robustly
accurate.
Figure 10. SE(t) for our standard twin experiment described in detail earlier when we utilize drifter information, and when we do not utilize drifter information. When the number of observations of height is L = 208, we see that without drifter information (blue line)
there is no synchronization and correspondingly inaccurate predictions (not shown). When information from 20 Lagrangian drifters is added during data assimilation using time delay nudging, SE(t) decreases very rapidly (red line), indicating predictions will be very accurate (also not shown). The efcacy of small numbers of drifters is clear in this example.
Figure 9. Initial positions for left panel ND = 20 drifters and right
panel ND = 64 drifters.
data remains to be seen however. Furthermore, the time delay method provides a natural way to incorporate this information into the analysis.
6 Discussion and summary
The transfer of information from measurements of a chaotic dynamical system to a quantitative model of the system is impeded when the number of measurements at each measurement time is below an approximate threshold Ls, which can be established in a twin experiment. Whartenby et al. (2013) previously showed that for a nonlinear model for shallow water ow, a standard nudging technique given by Eq. (2) requires direct observation of roughly 70 % of the dynamical
variables [notdef]h(r,t),u(r,t),v(r,t)[notdef] at each measurement time
to synchronize the model output with the observations.
Here we have demonstrated how information in the time delays of the observations may be used to reduce this requirement to about 30 %, in which only the height elds need be observed. Moreover, it appears Ls can be even further reduced by adding positional information from drifters, which interpolate the height eld at locations between grid points.
Although all this has been done on a simplied model of shallow water ow, implemented with only D = 3N2[Delta1] = 768
degrees of freedom, the process can be used to analyze increasingly realistic and complex models of coupled earth systems. Since the successful analysis of simulated data is a prerequisite for success with real data, this methodology provides a way to assess where one stands with respect to critical observability limits of the system at hand.
Furthermore, we expect that this formalism will generalize to systems substantially larger than the one presented here, although we do not underestimate the numerical challenges involved in its extension to, say, the scale of operational NWP models. We also suspect this issue of insufcient measurements to be a critical limitation in our current ability to predict the behavior of complex, chaotic systems. Since such systems are quite typical in practice, these issues need to be examined with more realistic models.
Harking back to the introduction, we note that the report by Cardinali (2013) indicates that 3040 million daily ob-
Nonlin. Processes Geophys., 24, 922, 2017 www.nonlin-processes-geophys.net/24/9/2017/
Z. An et al.: Estimating the state of a geophysical system with sparse observations 21
References
Abarbanel, H. D. I.: Analysis of Observed Chaotic Data, Springer,New York, 1996.
Abarbanel, H. D. I.: Predicting the Future: Completing Models of Observed Complex Systems, Springer-Verlag, New York, 2013.
Abarbanel, H. D. I., Creveling, D. R., Farsian, R., and Kostuk, M.: Dynamical State and Parameter Estimation, SIAM J. Appl. Dyn.
Syst., 8, 13411381, 2009.
Aeyels, D.: Generic observability of differentiable systems, SIAMJ. Control Optim., 19, 595603, 1981a.
Aeyels, D.: On the number of samples necessary to achieve observ-ability, Syst. Control Lett., 1, 9294, 1981b.
Blum, J.: Data assimilation for geophysical problems: variational and sequential techniques, Online Presentation from University of Nice Sophia Antipolis, France, 2010.
Blum, J., Le Dimet, F.-X., and Navon, I. M.: Chapter in Computational Methods for the Atmosphere and the Oceans, Volume 14: Special Volume of Handbook of Numerical Analysis, edited by: Temam, R. and Tribbia, J., Elsevier Science Ltd, New York, Hardbound, 784 pp., 2009.
Cardinali, C.: Data Assimilation: Observation Impact on the Short
Range Forecast, ECMWF Lecture Notes, available at: http://www.ecmwf.int/publications/
Web End =http:// http://www.ecmwf.int/publications/
Web End =www.ecmwf.int/publications/ , last access: January 2017, 2013. Courtier, P., Thpaut, J. N., and Hollingsworth, A.: A strategy for operational implementation of 4D-Var, using an incremental approach, Q. J. Roy. Meteor. Soc., 120, 13671387, 1994.
www.nonlin-processes-geophys.net/24/9/2017/ Nonlin. Processes Geophys., 24, 922, 2017
Figure 11. Comparison of the estimated and predicted elds [notdef]h(t),u(t),v(t)[notdef] between the truth (left column) and analyses, run with obser
vations of 128 height variables, both with (center column) and without (right column) drifters. Snapshots are taken 3 min (upper row) into the assimilation window, at 30 min at the end of the assimilation window (center row), and 90 min into the prediction window (bottom row).
servations are now available at the ECMWF, and that many
NWP models comprise upwards of 108 degrees of freedom. If the qualitative trends shown here, in which time delays provide successful predictions with only 30 % of the state variables observed, can be extended to substantially larger systems, then this method may indeed be useful for improving the forecasts of existing operational NWP models.
7 Data availability
Data used in this article are available at public data repository https://github.com/JasonAn/DataforNPG-23-1-2016
Web End =https://github.com/JasonAn/DataforNPG-23-1-2016 https://github.com/JasonAn/DataforNPG-23-1-2016
Web End = .
Acknowledgement. This work was funded in part under a grant from the US National Science Foundation (PHY-0961153). Partial support from the Department of Energy CSGF program (DE-FG02-97ER25308) for Daniel Rey is appreciated. Partial support from the MURI program (N00014-13-1-0205) sponsored by the Ofce of Naval Research is also acknowledged. We would also like to thank the reviewers for their thorough reading and thoughtful suggestions that have helped signicantly improve this manuscript.
Edited by: Z. Toth
Reviewed by: S. G. Penny and one anonymous referee
22 Z. An et al.: Estimating the state of a geophysical system with sparse observations
Durand, F., Gourdeau, L., Delcroix, T., and Verron, J.: Assimilation of sea surface salinity in a tropical Oceanic General Circulation Model (OGCM): A twin experiment approach, J. Geophys. Res., 107, 8004, doi:http://dx.doi.org/10.1029/2001JC000849
Web End =10.1029/2001JC000849 http://dx.doi.org/10.1029/2001JC000849
Web End = , 2002.
ECMWF: Lectures of Lars Isaksen and Stephen English at the European Centre for Medium-Range Weather Forecasts Training Course, June 2013.
Kantz, H. and Schreiber, T.: Nonlinear Time Series Analysis, 2ndEdn., Cambridge University Press, Cambridge, UK, 2004.
Kuznetsov, L., Ide, K., and Jones, C. K. R. T.: A Method for Assimilation of Lagrangian Data, Mon. Weather Rev., 131, 22472260, 2003.
Ma, R.: On the dimension of the compact invariant sets of certain nonlinear maps, in: Dynamical Systems and Turbulence, edited by: Rand, D. A. and Young, L.-S., Springer-Verlag, Lecture Notes in Mathematics, 898, 230242, 1981.
Mariano, A. J., Griffa, A., Zgkmen, T. M., and Zambianchi, E.: Lagrangian Analysis and Predictability of Coastal and Ocean Dynamics, J. Atmos. Ocean. Tech., 19, 11141126, 2002.Molcard, A., Piterbarg, L. I., Griffa, A., zgkmen, T. M., and Mariano, A. J.: Assimilation of drifter observations for the reconstruction of the Eulerian circulation eld, J. Geophys. Res., 108, 3056, doi:http://dx.doi.org/10.1029/2001JC001240
Web End =10.1029/2001JC001240 http://dx.doi.org/10.1029/2001JC001240
Web End = , 2003.
Lewis, J. M. and Derber, J. C.: The use of adjoint equations to solve a variational adjustment problem with advective constraints, Tellus A, 37A, 309322, 1985.
Lorenz, E. N.: Deterministic Nonperiodic Flow, J. Atmos. Sci., 20,130141, 1963.
Pazo, D., Carrassi, A., and Lopez, J. M.: Data-assimilation by delay-coordinate nudging, Q. J. Roy. Meteor. Soc., 142, 12901299, 2016.
Palatella, L., Carrassi, A., and Trevisan, A.: Lyapunov vectors and assimilation in the unstable subspace: theory and applications, J. Phys. A-Math. Theor., 46, 250301, doi:http://dx.doi.org/10.1088/1751-8113/46/25/250301
Web End =10.1088/1751- http://dx.doi.org/10.1088/1751-8113/46/25/250301
Web End =8113/46/25/250301 , 2013.
Pecora, L. and Carroll, T.: Synchronization in chaotic systems,Phys. Rev. Lett., 64, 821824, 1990.
Pedlosky, J.: Geophysical Fluid Dynamics, 2nd Edn., Springer-Verlag, 1987.
Pires, C., Vautard, R., and Talagrand, O.: On extending the limits of variational assimilation in nonlinear chaotic systems, Tellus A, 48, 96121, 1996.
Piterbarg, L. I.: Optimal estimation of Eulerian velocity eld given Lagrangian observations, Appl. Math. Model., 32, 21332148, 2008.
Press, W. H., Teukolsky, S. A., Vetterling, W. T., and Flannery, B. P.: Numerical Recipes, The Art of Scientic Computing, 3rd Edn., Cambridge University Press, 136139, 2007.
Quinn, J. C. and Abarbanel, H. D.: State and parameter estimation using Monte Carlo evaluation of path integrals, Q. J. Roy. Meteor. Soc., 136, 18551867, 2010.
Rabier, F., Jrvinen, H., Klinker, E., Mahfouf, J. F., and Simmons, A.: The ECMWF operational implementation of four-dimensional variational assimilation. I: Experimental results with simplied physics, Q. J. Roy. Meteor. Soc., 126, 11431170, 2000.
Rey, D., Eldridge, M., Kostuk, M., Abarbanel, H. D. I., Schumann-Bischoff, J., and Parlitz, U.: Accurate State and Parameter Estimation in Nonlinear Systems with Sparse Observations, Phys. Lett. A, 378, 869873, 2014a.
Rey, D., Eldridge, M., Morone, U., Abarbanel, H. D. I., Parlitz, U., and Schumann-Bischoff, J.: Using waveform information in nonlinear data assimilation, Phys. Rev. E, 90, 062916, doi:http://dx.doi.org/10.1103/PhysRevE.90.062916
Web End =10.1103/PhysRevE.90.062916 http://dx.doi.org/10.1103/PhysRevE.90.062916
Web End = , 2014b.
Sadourny, R.: The Dynamics of Finite Difference Models of the
Shallow Water Equations, J. Atmos. Sci., 32, 680689, 1975. Salman, H., Kuznetsov, L., Jones, C. K. R. T., and Ide, K.: A Method for Assimilating Lagrangian Data into a Shallow-Water-Equation Ocean Model, Mon. Weather Rev., 134, 10811101, 2006. Sauer, T., Yorke, J. A., and Casdagli, M.: Embedology, J. Stat.
Phys., 65, 579616, 1991.
Takens, F.: Detecting strange attractors in turbulence, in: Dynamical Systems and Turbulence, edited by: Rand, D. A. and Young, L.-S., Springer-Verlag, Lecture Notes in Mathematics, 898, 366 381, 1981.
Talagrand, O. and Courtier, P.: Variational Assimilation of Meteorological Observations With the Adjoint Vorticity Equation. I: Theory, Q. J. Roy. Meteor. Soc., 113, 13111328, 1987.
Thomson, R. E. and Emery, W. J.: Statistical Methods and Error Handling, Chapter 3, in: Data Analysis Methods in Physical Oceanography, 3rd Edn., Elsevier B. V., 219311, 2014. Trevisan, A., DIsidoro, M., and Talagrand, O.: Four dimensional variational assimilation in the unstable subspace and the optimal subspace dimension, Q. J. Roy. Meteor. Soc., 136, 487496, 2010.
Whartenby, W., Quinn, J., and Abarbanel, H. D. I.: The Number of Required Observations in Data Assimilation for a Shallow Water Flow, Mon. Weather Rev., 141, 25022518, 2013.
Nonlin. Processes Geophys., 24, 922, 2017 www.nonlin-processes-geophys.net/24/9/2017/
You have requested "on-the-fly" machine translation of selected content from our databases. This functionality is provided solely for your convenience and is in no way intended to replace human translation. Show full disclaimer
Neither ProQuest nor its licensors make any representations or warranties with respect to the translations. The translations are automatically generated "AS IS" and "AS AVAILABLE" and are not retained in our systems. PROQUEST AND ITS LICENSORS SPECIFICALLY DISCLAIM ANY AND ALL EXPRESS OR IMPLIED WARRANTIES, INCLUDING WITHOUT LIMITATION, ANY WARRANTIES FOR AVAILABILITY, ACCURACY, TIMELINESS, COMPLETENESS, NON-INFRINGMENT, MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. Your use of the translations is subject to all use restrictions contained in your Electronic Products License Agreement and by using the translation functionality you agree to forgo any and all claims against ProQuest or its licensors for your use of the translation functionality and any output derived there from. Hide full disclaimer
Copyright Copernicus GmbH 2017
Abstract
The problem of forecasting the behavior of a complex dynamical system through analysis of observational time-series data becomes difficult when the system expresses chaotic behavior and the measurements are sparse, in both space and/or time. Despite the fact that this situation is quite typical across many fields, including numerical weather prediction, the issue of whether the available observations are "sufficient" for generating successful forecasts is still not well understood. An analysis by Whartenby et al. (2013) found that in the context of the nonlinear shallow water equations on a β plane, standard nudging techniques require observing approximately 70% of the full set of state variables. Here we examine the same system using a method introduced by Rey et al. (2014a), which generalizes standard nudging methods to utilize time delayed measurements. We show that in certain circumstances, it provides a sizable reduction in the number of observations required to construct accurate estimates and high-quality predictions. In particular, we find that this estimate of 70% can be reduced to about 33% using time delays, and even further if Lagrangian drifter locations are also used as measurements.
You have requested "on-the-fly" machine translation of selected content from our databases. This functionality is provided solely for your convenience and is in no way intended to replace human translation. Show full disclaimer
Neither ProQuest nor its licensors make any representations or warranties with respect to the translations. The translations are automatically generated "AS IS" and "AS AVAILABLE" and are not retained in our systems. PROQUEST AND ITS LICENSORS SPECIFICALLY DISCLAIM ANY AND ALL EXPRESS OR IMPLIED WARRANTIES, INCLUDING WITHOUT LIMITATION, ANY WARRANTIES FOR AVAILABILITY, ACCURACY, TIMELINESS, COMPLETENESS, NON-INFRINGMENT, MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. Your use of the translations is subject to all use restrictions contained in your Electronic Products License Agreement and by using the translation functionality you agree to forgo any and all claims against ProQuest or its licensors for your use of the translation functionality and any output derived there from. Hide full disclaimer