<figure><img src="image_54.png" /><figcaption>FIG. 4. Lorenz-96 trajectories obtained using analog data assimilation procedures with the locally linear forecasting strategy, when only 20 variables are observed every 0.20 time steps. (top left) True simulation of the model with 40 variables, (top right) noisy and partial observations, (bottom left) reconstructed state trajectories via the AnEnKS with global analogs, and (bottom right) reconstructed state trajectories via the AnEnKS with local analogs [taking into account the 5 ($\nu$ = 2) nearest state components]. Only 10 Lorenz-96 cycles are shown for better visibility.</figcaption></figure>

for these three set of parameters, with $10^3$ Lorenz time steps each. Merging these three catalogs into a global catalog, we apply the proposed AnDA using as observations the true integration resulting from Lorenz-63 model with $\theta_1$ parameter values. As a by-product of the analog strategy, we can infer the underlying model parameterization from the observed partial observations. The reported experiments (Fig. 6) apply the AnPF procedure with the locally constant analog method and a multinomial sampling scheme using 100 particles. Such a choice was motivated by the desire of keeping track of the particles and their source catalog, which is harder to achieve with the other AnDA algorithms, since the particles would be elements from the catalog and the AnPF assigns a weight to each particle. This makes it easier to select at each time the particle with the biggest weight and to know from which catalog it came from.

At every assimilation time step, we determine which parameterization most ensemble members come from,

and then calculate the proportion of the presence of each parameterization. As expected, the true parameterization (red, parameterization $\theta_1$) is more represented. The proportions for $\theta_1$, $\theta_2$, and $\theta_3$ are around,

<table><caption>TABLE 1. RMSE of the reconstruction of Lorenz-96 state evolution using different forecasting strategies and data assimilation techniques. The catalog size corresponds to $10^3$ Lorenz-96 times (equivalent to 13 yr) and the number of members/particles is $N = 1000$.</caption><thead><tr><th>Method</th><th>Locally constant</th><th>Locally incremental</th><th>Locally linear</th></tr></thead><tbody><tr><td>Gaussian</td><td></td><td></td><td></td></tr><tr><td>AnEnKF</td><td>1.826</td><td>1.785</td><td>1.403</td></tr><tr><td>AnPF</td><td>3.174</td><td>4.224</td><td>4.4616</td></tr><tr><td>AnEnKS</td><td>1.320</td><td>1.287</td><td>0.970</td></tr><tr><td>Multinomial</td><td></td><td></td><td></td></tr><tr><td>AnEnKF</td><td>1.814</td><td>1.774</td><td>1.413</td></tr><tr><td>AnPF</td><td>2.989</td><td>4.412</td><td>4.729</td></tr><tr><td>AnEnKS</td><td>1.313</td><td>1.288</td><td>1.093</td></tr></tbody></table>
<figure><img src="image_55.png" /><figcaption>FIG. 5. Reconstruction of Lorenz-63 trajectories for different catalog sizes in the analog data assimilation pro- cedures, when only the first component of the state is observed every 0.08 time steps. (left) RMSE as a function of the size of the catalog for different analog data assimilation strategies: AnEnKF (green), AnPF (blue), and AnEnKS (red). For benchmarking purposes, data assimilation results with true Lorenz-63 equations are given in straight lines. (right) Time series of the first component of the true state (black solid line), associated noisy ob- servations (black asterisks), mean reconstructed series (solid lines), and 10 analyzed members/particles (dashed lines) with analog data assimilation strategies, namely AnEnKF (green), AnPF (blue), and AnEnKS (red), using a catalog of 10³ Lorenz-63 times (equivalent to 8 yr).</figcaption></figure>

60%, 16%, and 24%, respectively, proving the ability of the methodology to detect the source of the noisy and partial observation (here, only coming from $\theta_1$). To analyze the results more thoroughly, we calculate the RMSE of the reconstruction using (i) the three catalogs as shown before, (ii) only the good catalog, and (iii) only the two “bad” catalogs. The RMSEs are (i) 1.287, (ii) 1.207, and (iii) 1.424, respectively. These results show that having other catalogs with

different parameterization degrade the RMSE but the filter is still performing well. This experiment gives insights on the problem of the assimilation of variables that may switch between different dynami- cal modes. Analog data assimilation can deal with this problem in a simpler manner than classical data as- similation, through the concatenation of the catalogs issued from different parameterizations into a single catalog.

<figure><img src="image_56.png" /><figcaption>FIG. 6. Identification of Lorenz-63 model parameterizations using a multiparameterization catalog in the analog data assimilation, when only the first component of the state is observed every 0.08 time step. (left) Examples of Lorenz-63 trajectories generated with three different parameterizations: $\theta_1 = (10, 28, 8/3)$ (red), $\theta_2 = (7, 28, 8/3)$ (blue), and $\theta_3 = (13, 28, 8/3)$ (green). (right) Result of the AnPF on the first Lorenz-63 variable using the three catalogs associated with parameterizations $\{\theta_i\}_{1,2,3}$ for $3 \times 10^3$ Lorenz-63 times (equivalent to $3 \times 8$ yr) when only observations from parameterization $\theta_1 = (10, 28, 8/3)$ are provided. The figure shows the AnPF particles trajectories (blue), the AnPF result (red), and the true trajectory (green).</figcaption></figure>
<figure><img src="image_57.png" /><figcaption>FIG. 7. Results of the reconstruction of Lorenz-63 trajectories from noisy catalogs: (left) examples of noisy Lorenz-63 trajectories for different noise levels: $\psi_1^2 = 0.5$ (red), $\psi_2^2 = 1$ (blue), and $\psi_3^2 = 2$ (green). (right) Results of the AnEnKS using noisy catalogs corresponding to $10^3$ Lorenz-63 times (equivalent to 8 yr) when only observations with variance $R = 2$ are provided. We also plot the 95% confidence interval computed from the smoothing covariances.</figcaption></figure>

### 3) EXPERIMENT 3

Whereas previous experiments consider catalogs produced from noise-free trajectories, here we evaluate the sensitivity of the AnDA procedures when the catalog may involve noisy trajectories of the considered system. Acquisition systems typically involve such noise patterns, which may relate for instance to both environmental constraints and measurement uncertainties. We simulate noisy catalogs for Lorenz-63 dynamics as follows: we artificially degrade the transition between consecutive states with a Gaussian additive noise. We performed experiments with different noise variances $\psi^2 = \{0.5, 1, 2\}$ to evaluate the sensitivity of AnDA procedures with respect to the signal-to-noise ratio. As illustrated in Fig. 7, the trajectories of these experiments are extremely noisy. Table 2 reports the RMSE of the different AnDA algorithms with the locally linear analog forecasting operator and 100 ensemble members/particles. As expected, the RMSE increases with the variance of the additive noise. The AnEnKS clearly outperforms the other AnDA algorithms, which highlights its greater robustness. Figure 7 further illustrates that the AnEnKS is able to correctly track the true state of the system, even for highly degraded catalogs ($\psi^2 = 2$, green curve). For a high signal-to-noise ratio (i.e., low perturbations) ($\psi^2 = 0.5$, red curve), reconstructed trajectories are very close to the ones obtained with a noise-free catalog.

### 6. Conclusions and perspectives

The present paper demonstrates the potential of data-driven schemes for data assimilation. We propose and evaluate efficient yet simple data-driven forecasting

strategies that can be coupled with classical stochastic filters (viz., the ensemble Kalman filter/smoothen and the particle filter). We set a unified framework that we call the analog data assimilation (AnDA). The key features of the AnDA are twofold: (i) it relies on a data-driven representation of the state dynamics, and (ii) it does not require online evaluations of dynamical models based on physical equations. The relevance of the AnDA is tangible when the dynamical system of interest demands tremendous and time-consuming physical modeling efforts and/or uncertainties are difficult to assess. In cases when large observational or model-simulated datasets of the considered system are available, AnDA can both support or compete with classical data assimilation schemes. As a proof concept, we demonstrate the relevance of the proposed methodology to retrieve the chaotic behavior of the Lorenz-63 and Lorenz-96 models. We performed numerical experiments to evaluate critical aspects of the method, especially the relevant combinations of analog forecasting strategies and of stochastic filters as well as the exploitation of noisy and noise-free catalogs.

<table><caption>TABLE 2. RMSE of the reconstruction of Lorenz-63 trajectories from noisy catalogs: we vary the variance of an additive Gaussian noise in the creation of the catalogs and apply analog data assimilation procedures with the locally linear operator with a catalog size of $10^3$ Lorenz-63 times, when only the first component of the state is observed every 0.08 time step with observation noise variance $R = 2$.</caption><thead><tr><th>Method</th><th>ψ<sub>1</sub><sup>2</sup> = 0.5</th><th>ψ<sub>2</sub><sup>2</sup> = 1</th><th>ψ<sub>3</sub><sup>2</sup> = 2</th></tr></thead><tbody><tr><td>AnEnKF</td><td>1.926</td><td>2.136</td><td>2.681</td></tr><tr><td>AnPF</td><td>1.652</td><td>1.961</td><td>2.313</td></tr><tr><td>AnEnKS</td><td>1.233</td><td>1.561</td><td>2.142</td></tr></tbody></table>
All the reported experiments were carried out using the AnDA Python/Matlab library (https://github.com/ptandeo/AnDA), which includes the Lorenz-63 and Lorenz-96 systems. In the spirit of reproducible research, the user can conduct the different experiments shown in this paper.

Overall, the reported results demonstrate the relevance of the proposed analog data assimilation methods, even with highly damaged catalogs. They suggest that AnEnKS combined with locally incremental or locally linear analog forecasting leads to the best reconstruction performance, the locally incremental version being the most robust to noisy settings. Moreover, the flexibility of the analog data assimilation demonstrates the potential for the identification of hidden underlying dynamics from a series of partial observations.

The main pillar of our data-driven approach is the catalog. As such, analog data assimilation deeply relates to the quality and representativity of the catalog. In our experiments, we assumed that we were provided with large-scale catalogs of complete states of the system of interest. While catalogs built from numerical simulations fulfill this assumption, observational datasets (e.g., satellite remote sensing or in situ data) typically involve missing data, which may require specific strategies to be dealt with in the building of the catalogs. In this respect, local analogs obviously appear much more flexible than global ones, as partial observations provide relevant exemplars for the creation of catalogs for local analogs.

The application of analog data assimilation to high-dimensional systems is another future challenge. As detailed in Van den Dool (1994), the number of elements in a catalog shall grow exponentially with the intrinsic dimension of the state to guarantee the retrieval of analogs at a given precision. This makes unrealistic the direct application of analog strategies to state space with an intrinsic dimensionality above 10. As a consequence, global analog forecasting operators are most likely inappropriate for high-dimensional systems. By contrast, local analogs provide a means to decompose the analog forecasting of the high-dimensional state into a series of local and low-dimensional analog forecasting operations. This is regarded as the key explanation for the much better performance reported for the local analog data assimilation for Lorenz-96 dynamics using catalogs of about a million of exemplars (Fig. 4). For real-world applications to high-dimensional systems, for instance to ocean and atmosphere dynamics, the combination of such local analog strategies to multiscale decompositions (Mallat 1989) arise as a promising research direction as illustrated in Fablet et al. (2017). Such multiscale decompositions are expected to enhance the spatial redundancy, with a view to

building the requested catalogs of millions to hundreds of millions of exemplars (for an intrinsic dimensionality between 4 and 7, see the appendix) from observation or simulation datasets over a few decades. Another important aspect that controls the effective size of the catalog is the evolution of the system in time. The more nonlinear the dynamics, the greater the number of requested exemplars in the global catalog to learn the forecast operator and the spread of the prediction.

We believe that this study opens new research avenues for the analysis, reconstruction, and understanding of the dynamics of geophysical systems using data-driven techniques. Such techniques will benefit from the increasing availability of large-scale historical observational and/or simulated datasets. Beyond the wide range of possible applications, future research should further investigate methodological issues. First of all, our study demonstrates the relevance of the analog particle filter, but as mentioned in section 5, the AnPF suffers from degeneracy and sample impoverishment. We may point out that complementary experiments with particle smoother schemes (not shown in this paper) resulted in numerical instabilities. The derivation of the analog particle smoother then remains an open question. In addition to advanced particle filters as proposed in Van Leeuwen (2010) and Pitt and Shephard (1999), one might also benefit from the straightforward applications of the analog procedure in reverse time, which is not generally possible for model-driven schemes. A second direction for future work lies in the design of the kernel used by the analog forecasting operators. Whereas we considered a Gaussian kernel, other kernels have been proposed in the literature; for example, using Procrustes distance instead of the Euclidean distance (McDermott and Wikle 2016) or different weighing strategies (Delle Monache et al. 2011). The explicit derivation of the mapping associated with a kernel as considered in Zhao and Giannakis (2014) may also be a promising alternative to state the analog data assimilation in a kernel-derived lower-dimensional space. The theoretical characterization of the asymptotic behavior of analog data assimilation schemes is also an interesting avenue of research. Similarly to the theoretical analysis of ensemble Kalman filters and particle filters (Le Gland et al. 2009), the derivation of convergence conditions, possibly associated with reconstruction bounds, would be of key interest to bound the reconstruction performance of the proposed analog schemes with respect to their model-driven counterpart.

*Acknowledgments.* We thank all the researchers from various fields who provided careful and constructive comments on the original paper especially Bertrand
Chapron, Valérie Monbet, and Anne Cuzol. The authors would also like to thank Phi Viet Huynh for his valuable contribution to both the AnDA Matlab toolbox and the AnDA Python library. We are also grateful to the two anonymous reviewers, whose comments helped to improve the manuscript. We thank Geraint Jones for his English grammar corrections. This work was supported by ANR (Agence Nationale de la Recherche, Grant ANR-13-MONU-0014), Labex Cominlabs (Grant SEACS), the Brittany council, and a “Futur et Ruptures” postdoctoral grant from Institut Mines-Télécom.

# APPENDIX

## Operational Count of the AnDA Applied for High-Dimensional Applications

This appendix aims at giving an estimate of the operations involved when applying the AnDA for a realistic large-scale application. We discuss the computational cost of the analog forecasting, which is specific to the AnDA. The latter directly relates to the cost of the K-nearest neighbor (K-NN) step.

In case of large-scale catalogs, an exhaustive search strategy is not suitable and the use of space-partitioning data structures, the most popular ones being K-d trees (Bentley 1975) and Ball trees (Omohundro 1989), appears necessary. These structures speed up the K-NN search, at the expense of an approximate search for nearest neighbors. Let us denote by D the dimension of the system of interest. Making a choice between K-d trees or ball trees depends mostly on the dimensionality of the system. The K-d trees are known to perform well in dimensions D < 20, while ball trees are more suitable to dimensions higher than 20 but come with a high cost of space partitioning (Witten et al. 2016). In this appendix we focus on the use of K-d trees, which are natural candidates for local analogs with a small component-wise local neighborhood v or using a preliminary dimensionality reduction algorithm (such as empirical orthogonal functions). A comparison between K-d trees and ball trees is out of the scope of this work.

Let $N_{\text{data}}$ be the size of the catalog (the number of samples from where to look for analogs), and K the number of nearest neighbors to be retrieved. Let us recall that v is the size of the local neighborhood used for the search for local analogs. Van den Dool (1994) derived a relationship between the local neighborhood size and the amount of the data needed to find an analog with a given precision. With the assumption that the components of the states follow a multivariate Gaussian distribution and have the same variance $sd^2$, finding K samples that have a distance lower than $\epsilon$ for all

the components of the neighborhood with a probability of 95%, needs the number of data to be on average as follows:

* Global analogs:

$$N_{\text{global}} \geq K \frac{\ln(0.05)}{\ln(1 - \alpha^D)} \approx \frac{3K}{\alpha^D}, \quad (A1)$$

* Local analogs:

$$N_{\text{local}} \geq K \frac{\ln(0.05)}{\ln(1 - \alpha^{2\nu+1})} \approx \frac{3K}{\alpha^{2\nu+1}}, \quad (A2)$$

where $\alpha$ is the integral of the standard Gaussian probability density function from $-\epsilon/(\sqrt{2}sd)$ to $-\epsilon/(\sqrt{2}sd)$.

We present now the operational count for one ensemble member (or particle) involved in the forecasting, for both global and local analogs. In each case, we distinguish the computational cost of the creation of the K-d trees and the search of K nearest neighbors:

* Global analogs:

* Creation of the K-d tree: $O[DN_{\text{global}}\log(N_{\text{global}})]$.

* Search for K global analogs: $O[KD\log(N_{\text{global}})]$.

* Local analogs:

* Creation of D K-d trees (for every dimension in D): $O[D(2\nu + 1)N_{\text{local}}\log(N_{\text{local}})]$.

* Search for K local analogs of component-wise neighborhood v: $O[DK(2\nu + 1)\log(N_{\text{local}})]$.

Note that using local analogs requires constructing a K-d tree for every dimension in D. Construction of the K-d trees can be done offline (1 “big” K-d tree for the global strategy and D “small” K-d trees for the local strategy), then the cost of these construction can be amortized over the high number of queries that needs to be answered during analog data assimilation. However, in terms of memory storage, storing a global K-d tree could be prohibitive, contrarily to small local K-d trees that can be created, used, then freed for the creation of the next K-d tree of the next dimension (if there is no sufficient memory to stock D small local K-d trees). Keep in mind that we need to have $(2\nu + 1) \ll D$ for local analogs to be of relevance.

Let us take an example using the Lorenz-96 model: $D = 40$, $\nu = 2$. Looking for $K = 50$ analogs, with an $\alpha = 0.15$ we would need $N_{\text{global}} \approx 10^{35}$, which is very prohibitive; however, we would only need $N_{\text{local}} \approx 2 \times 10^6$ samples using local analogs.

# REFERENCES

Anderson, J. L., 2001: An ensemble adjustment Kalman filter for data assimilation. Mon. Wea. Rev., 129, 2884–2903, doi:10.1175/1520-0493(2001)129<2884:AEAKFF>2.0.CO;2.
—, 2007: Exploring the need for localization in ensemble data assimilation using a hierarchical ensemble filter. *Physica D*, **230**, 99–111, doi:10.1016/j.physd.2006.02.011.

—, 2012: Localization and sampling error correction in ensemble Kalman filter data assimilation. *Mon. Wea. Rev.*, **140**, 2359–2371, doi:10.1175/MWR-D-11-00013.1.

—, and S. L. Anderson, 1999: A Monte Carlo implementation of the nonlinear filtering problem to produce ensemble assimilations and forecasts. *Mon. Wea. Rev.*, **127**, 2741–2758, doi:10.1175/1520-0493(1999)127<2741:AMCIOT>2.0.CO;2.

Atencia, A., and I. Zawadzki, 2015: A comparison of two techniques for generating nowcasting ensembles. Part II: Analogs selection and comparison of techniques. *Mon. Wea. Rev.*, **143**, 2890–2908, doi:10.1175/MWR-D-14-00342.1.

Bentley, J. L., 1975: Multidimensional binary search trees used for associative searching. *Commun. ACM*, **18**, 509–517, doi:10.1145/361002.361007.

Bocquet, M., C. A. Pires, and L. Wu, 2010: Beyond Gaussian statistical modeling in geophysical data assimilation. *Mon. Wea. Rev.*, **138**, 2997–3023, doi:10.1175/2010MWR3164.1.

Bowler, N. E., J. Flowerdew, and S. R. Pring, 2013: Tests of different flavours of EnKF on a simple model. *Quart. J. Roy. Meteor. Soc.*, **139**, 1505–1519, doi:10.1002/qj.2055.

Burgers, G., P. Jan van Leeuwen, and G. Evensen, 1998: Analysis scheme in the ensemble Kalman filter. *Mon. Wea. Rev.*, **126**, 1719–1724, doi:10.1175/1520-0493(1998)126<1719:ASITEK>2.0.CO;2.

Chin, T., M. Turmon, J. Jewell, and M. Ghil, 2007: An ensemble-based smoother with retrospectively updated weights for highly nonlinear systems. *Mon. Wea. Rev.*, **135**, 186–202, doi:10.1175/MWR3353.1.

Cleveland, W. S., 1979: Robust locally weighted regression and smoothing scatterplots. *J. Amer. Stat. Assoc.*, **74**, 829–836, doi:10.1080/01621459.1979.10481038.

Cosme, E., J. Verron, P. Brasseur, J. Blum, and D. Auroux, 2012: Smoothing problems in a Bayesian framework and their linear Gaussian solutions. *Mon. Wea. Rev.*, **140**, 683–695, doi:10.1175/MWR-D-10-05025.1.

Delle Monache, L., T. Nipen, Y. Liu, G. Roux, and R. Stull, 2011: Kalman filter and analog schemes to postprocess numerical weather predictions. *Mon. Wea. Rev.*, **139**, 3554–3570, doi:10.1175/2011MWR3653.1.

—, I. Djalalova, and J. Wilczak, 2014: Analog-based postprocessing methods for air quality forecasting. *Air Pollution Modeling and Its Application XXIII*, D. Steyn and R. Mathur, Eds., Springer, 237–239, doi:10.1007/978-3-319-04379-1_38.

Dormand, J. R., and P. J. Prince, 1980: A family of embedded Runge-Kutta formulae. *J. Comput. Appl. Math.*, **6**, 19–26, doi:10.1016/0771-050X(80)90013-3.

Evensen, G., 2007: *Data Assimilation: The Ensemble Kalman Filter*. Springer-Verlag, 280 pp., doi:10.1007/978-3-540-38301-7.

—, and P. J. Van Leeuwen, 2000: An ensemble Kalman smoother for nonlinear dynamics. *Mon. Wea. Rev.*, **128**, 1852–1867, doi:10.1175/1520-0493(2000)128<1852:AEKSFN>2.0.CO;2.

Fablet, R., P. H. Viet, R. Lguensat, and B. Chapron, 2017: Data-driven assimilation of irregularly-sampled image time series. *IEEE Int. Conf. on Image Processing (ICIP 2017)*, Beijing, China, IEEE, WQ-PB.2.

Hamilton, F., T. Berry, and T. Sauer, 2016: Ensemble Kalman filtering without a model. *Phys. Rev. X*, **6**, 011021, doi:10.1103/PhysRevX.6.011021.

Hansen, B., 2000: Econometrics. Department of Economics, University of Wisconsin, 427 pp., http://www.ssc.wisc.edu/~bhansen/econometrics/Econometrics.pdf.

He-Guelton, L., R. Fablet, B. Chapron, and J. Tournadre, 2015: Learning-based emulation of sea surface wind fields from numerical model outputs and SAR data. *IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens.*, **8**, 4742–4750, doi:10.1109/JSTARS.2015.2496503.

Hong, S.-Y., and J. Dudhia, 2012: Next-generation numerical weather prediction: Bridging parameterization, explicit clouds, and large eddies. *Bull. Amer. Meteor. Soc.*, **93**, ES6-ES9, doi:10.1175/2011BAMS3224.1.

Hoteit, I., D.-T. Pham, G. Triantafyllou, and G. Korres, 2008: A new approximate solution of the optimal nonlinear filter for data assimilation in meteorology and oceanography. *Mon. Wea. Rev.*, **136**, 317–334, doi:10.1175/2007MWR1927.1.

—, X. Luo, and D.-T. Pham, 2012: Particle Kalman filtering: A nonlinear Bayesian framework for ensemble Kalman filters. *Mon. Wea. Rev.*, **140**, 528–542, doi:10.1175/2011MWR3640.1.

—, D.-T. Pham, M. Gharamti, and X. Luo, 2015: Mitigating observation perturbation sampling errors in the stochastic EnKF. *Mon. Wea. Rev.*, **143**, 2918–2936, doi:10.1175/MWR-D-14-00088.1.

Kalnay, E., 2003: *Atmospheric Modeling, Data Assimilation and Predictability*. Cambridge University Press, 345 pp.

Le Gland, F., V. Monbet, and V.-D. Tran, 2009: Large sample asymptotics for the ensemble Kalman filter. Research Rep. RR-7014, INRIA, 25 pp., https://hal.inria.fr/inria-00409060/document.

Lorenc, A., and Coauthors, 2000: The Met. Office global three-dimensional variational data assimilation scheme. *Quart. J. Roy. Meteor. Soc.*, **126**, 2991–3012, doi:10.1002/qj.49712657002.

Lorenz, E. N., 1969: Atmospheric predictability as revealed by naturally occurring analogues. *J. Atmos. Sci.*, **26**, 636–646, doi:10.1175/1520-0469(1969)26<636:APARBN>2.0.CO;2.

—, 1996: Predictability—A problem partly solved. *Proc. Seminar on Predictability*, Reading, United Kingdom, ECMWF, 18 pp., https://www.ecmwf.int/sites/default/files/elibrary/1995/10829-predictability-problem-partly-solved.pdf.

Mallat, S. G., 1989: A theory for multiresolution signal decomposition: The wavelet representation. *IEEE Trans. Pattern Anal. Mach. Intell.*, **11**, 674–693, doi:10.1109/34.192463.

McDermott, P. L., and C. K. Wikle, 2016: A model-based approach for analog spatio-temporal dynamic forecasting. *Environmetrics*, **27**, 70–82, doi:10.1002/env.2374.

Miller, R. N., M. Ghil, and F. Gauthiez, 1994: Advanced data assimilation in strongly nonlinear dynamical systems. *J. Atmos. Sci.*, **51**, 1037–1056, doi:10.1175/1520-0469(1994)051<1037:ADAISN>2.0.CO;2.

Nagarajan, B., L. Delle Monache, J. P. Hacker, D. L. Rife, K. Searight, J. C. Knievel, and T. N. Nipen, 2015: An evaluation of analog-based postprocessing methods across several variables and forecast models. *Wea. Forecasting*, **30**, 1623–1643, doi:10.1175/WAF-D-14-00081.1.

Omohundro, S. M., 1989: Five balltree construction algorithms. International Computer Science Institute, Berkeley, CA, 22 pp., http://ftp.icsi.berkeley.edu/ftp/pub/techreports/1989/tr-89-063.pdf.

Ott, E., and Coauthors, 2004: A local ensemble Kalman filter for atmospheric data assimilation. *Tellus*, **56A**, 415–428, doi:10.3402/tellusa.v56i5.14462.

Pham, D. T., 2001: Stochastic methods for sequential data assimilation in strongly nonlinear systems. *Mon. Wea. Rev.*, **129**, 1194–1207, doi:10.1175/1520-0493(2001)129<1194:SMFSDA>2.0.CO;2.
Pitt, M. K., and N. Shephard, 1999: Filtering via simulation: Auxiliary particle filters. *J. Amer. Stat. Assoc.*, **94**, 590–599, doi:10.1080/01621459.1999.10474153.

Ruiz, J. J., M. Pulido, and T. Miyoshi, 2013: Estimating model parameters with ensemble-based data assimilation: A review. *J. Meteor. Soc. Japan*, **91**, 79–99, doi:10.2151/jmsj.2013-201.

Schenk, F., and E. Zorita, 2012: Reconstruction of high resolution atmospheric fields for northern Europe using analog-upscaling. *Climate Past*, **8**, 1681–1703, doi:10.5194/cp-8-1681-2012.

Schölkopf, B., and A. J. Smola, 2001: *Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Bayard*. MIT Press, 648 pp.

Tandeo, P., and Coauthors, 2015a: Combining analog method and ensemble data assimilation: Application to the Lorenz-63 chaotic system. *Machine Learning and Data Mining Approaches to Climate Science*, V. Lakshmanan et al., Eds., Springer, 3–12, doi:10.1007/978-3-319-17220-0_1.

—, M. Pulido, and F. Lott, 2015b: Offline parameter estimation using EnKF and maximum likelihood error covariance estimates: Application to a subgrid-scale orography parametrization. *Quart. J. Roy. Meteor. Soc.*, **141**, 383–395, doi:10.1002/qj.2357.

Van den Dool, H., 1994: Searching for analogues, how long must we wait? *Tellus*, **46A**, 314–324, doi:10.3402/tellusa.v46i3.15481.

Van Leeuwen, P. J., 2009: Particle filtering in geophysical systems. *Mon. Wea. Rev.*, **137**, 4089–4114, doi:10.1175/2009MWR2835.1.

—, 2010: Nonlinear data assimilation in geosciences: An extremely efficient particle filter. *Quart. J. Roy. Meteor. Soc.*, **136**, 1991–1999, doi:10.1002/qj.699.

Whitaker, J. S., and T. M. Hamill, 2002: Ensemble data assimilation without perturbed observations. *Mon. Wea. Rev.*, **130**, 1913–1924, doi:10.1175/1520-0493(2002)130<1913:EDAWPO>2.0.CO;2.

Wilby, R. L., and T. Wigley, 1997: Downscaling general circulation model output: A review of methods and limitations. *Prog. Phys. Geogr.*, **21**, 530–548, doi:10.1177/030913339702100403.

Witten, I. H., E. Frank, M. A. Hall, and C. J. Pal, 2016: *Data Mining: Practical Machine Learning Tools and Techniques*. 4th ed. Morgan Kaufmann, 654 pp.

Yiou, P., 2014: AnaWEGE: A weather generator based on analogues of atmospheric circulation. *Geosci. Model Dev.*, **7**, 531–543, doi:10.5194/gmd-7-531-2014.

Zhao, Z., and D. Giannakis, 2014: Analog forecasting with dynamics-adapted kernels. *Nonlinearity*, **29**, 2888–2939, doi:10.1088/0951-7715/29/9/2888.
## 4.4 Tandeo, Pulido et Lott (2015) [QJRMS]

**Contexte** Cet article correspond à mon principal travail de recherche lorsque j’étais postdoctorant en Argentine. Grâce à une collaboration entre le LMD et CONICET (agence nationale de recherche argentine), nous avons travaillé sur l’estimation de paramètres dans les processus sous-maille, modélisant des processus fortement non linéaires induits par l’orographie. Ce travail m’a permis d’attaquer un problème de modélisation physique complexe avec des outils d’assimilation de données. Ce travail a été publié dans QJRMS, un journal en sciences atmosphériques, avec une forte composante en assimilation de données. Les développements méthodologiques que j’ai proposés lors de cette étude sont toujours d’actualité et, avec mes collègues argentins, nous continuons à améliorer les méthodes d’estimation de paramètres.

**Résumé** Des travaux récents ont montré que les paramètres contrôlant les paramétrisations de processus physiques dans les modèles climatiques peuvent être estimés à partir d'observations utilisant des techniques de filtrage. Dans cet article, nous proposons une approche d'estimation des paramètres hors ligne, sans estimer l'état du modèle climatique. Elle est basée sur le filtre de Kalman d'ensemble et une estimation itérative des matrices de covariance d'erreur de l'ébauche et des observations, à l'aide d'un algorithme maximisant la vraisemblance. La technique est mise en œuvre dans un processus sous-maille orographique 1D vertical. Tout d'abord, la méthode d'estimation est évaluée à l'aide d'expériences ju-melles. Ensuite, la méthode est utilisée avec des observations synthétiques pour évaluer comment les paramètres du processus sous-maille changent lorsque la résolution orographique d'un modèle de circu-lation générale est augmentée. Notre analyse révèle que, lorsque la résolution orographique augmente, la modélisation sous-maille doit tenir compte de l'effet d'abri qui peut se produire à faible altitude entre des sommets montagneux situés dans un même point de grille.
# Offline parameter estimation using EnKF and maximum likelihood error covariance estimates: Application to a subgrid-scale orography parametrization

P. Tandeo,<sup>a,b*</sup> M. Pulido<sup>a</sup> and F. Lott<sup>c</sup>

Recent work has shown that the parameters controlling parametrizations of the physical processes in climate models can be estimated from observations using filtering techniques. In this article, we propose an offline parameter estimation approach, without estimating the state of the climate model. It is based on the Ensemble Kalman Filter (EnKF) and an iterative estimation of the error covariance matrices and of the background state using a maximum likelihood algorithm. The technique is implemented in a subgrid-scale orography (SSO) parametrization scheme which works in a single vertical column. First, the parameter estimation technique is evaluated using twin experiments. Then, the technique is used with synthetic observations to estimate how the parameters of the SSO scheme should change when the resolution of the input orography dataset of a general circulation model is increased. Our analysis reveals that, when the resolution of the orography dataset increases, the scheme should take into account the dynamical sheltering that can occur at low levels between mountain peaks located within the same gridbox area.

**Key Words:** offline parameter estimation; EnKF; EM algorithm; subgrid-scale orography parametrization

Received 21 June 2013; Revised 4 February 2014; Accepted 27 February 2014; Published online in Wiley Online Library

## 1. Introduction

Numerical models, including atmospheric/oceanic general circulation models (GCMs) and current earth system models, contain several physical parametrizations with a large number of parameters. Climate predictions using these numerical models are sensitive to the large set of parameters that are present in the physical parametrizations (cf. Stainforth et al., 2005). Most of these unknown physical parameters can not be determined directly from observations and are generally manually tuned. This subjective approach is excessively time demanding and does not give optimal results. Moreover, if the horizontal resolution of the model or of an input dataset is increased or a parametrization scheme is changed, the physical parameters need to be re-evaluated. To address these issues, several authors (e.g. Jackson et al., 2004; Severijns and Hazeleger, 2005) propose estimating the physical parameters objectively, defining a cost function based on the root mean square error (RMSE) criterion. The idea is to find the optimal set of parameter values that gives the minimum RMSE and produces the lowest model error. However, nonlinear model responses may produce multiple local minima in the cost function (cf. Posselt and Bishop, 2012), and thus sophisticated optimization algorithms are required to find the global minimum corresponding to the optimal parameters. Such optimization

algorithms are usually too expensive computationally to be employed in sophisticated models. An alternative consists of supposing that the parameters are stochastic and of estimating them using filtering techniques (e.g. Annan and Hargreaves, 2004; Posselt and Bishop, 2012; Ruiz et al., 2013; Schirber et al., 2013). The basic idea to estimate the parameters is based on an augmented state composed of both the state of the system and the physical parameters in a nonlinear Gaussian state-space model. This online estimation is a tough problem in practice. Even a simple linear state equation with multiplicative parameters behaves nonlinearly for parameter estimation (Yang and Delsole, 2009).

Another approach consists of estimating the physical parameters independently of the state of the system. The particular advantage of using an offline estimation technique is that the control space is reduced from $10^7$ to just a few dimensions. This drastic reduction in size permits us to conduct several model/parametrization evaluations as is often needed in parameter estimation. One disadvantage of offline techniques is that they cannot take into account the feedback of the changes that the parametrization produces onto the parametrization itself. Nevertheless, for a subgrid-scale orography (SSO) scheme, this issue should not be too critical, since most of the flow changes produced are advected downstream (for instance in the form of potential vorticity banners; Figure 13(c) in Lott, 1995).

<footnote><sup>a</sup>Department of Physics, Universidad Nacional del Nordeste, Corrientes, Argentina</footnote>

<footnote><sup>b</sup>Lab-STICC - Pôle CID, Telecom Bretagne, Brest, France</footnote>

<footnote><sup>c</sup>Laboratoire de Meteorologie Dynamique, Ecole Normale Superieure, Paris, France</footnote>

<footnote>*Correspondence to: P. Tandeo, Telecom Bretagne, 29280 Plouzané, France.
E-mail: pierre.tandeo@telecom-bretagne.eu</footnote>
Accordingly, the feedback can be neglected if the mountains considered are not close to the lee of other mountains.

In order to conduct an offline estimation of physical parameters, the parametrization should be compared to observations, for instance the Pyrenees Experiment (PYREX) campaign, in which surface drag and momentum fluxes were measured over a transect of the Pyrenees mountains (cf. Bougeault et al., 1990). In this case, the mountain massif can be considered to be entirely located within a model gridbox area of a climate model, so drag and momentum flux can be directly compared to the same quantities predicted by the scheme over the same area. Therefore, we can validate SSO schemes using a single vertical column. This approach is often used prior to the implementation of the schemes in GCMs (e.g. the offline tests of the scheme in single vertical columns using the PYREX data in Lott and Miller, 1997).

Pulido and Thuburn (2005, 2008) showed that a four-dimensional variational data assimilation technique can be used to estimate the missing momentum forcing due to the unresolved/subgrid-scale gravity waves in the stratosphere. This missing momentum forcing was used to estimate optimal parameters of a non-orographic gravity wave parametrization in Pulido et al. (2012). Using twin experiments, they showed that the variational data assimilation technique does not converge towards the optimal parameters because of the nonlinear response of the parametrization to parameter perturbations. They employed a time-demanding genetic algorithm to overcome these difficulties. In the present work, we propose a similar offline parameter estimation procedure, but using an ensemble-based data assimilation technique to estimate the optimal parameters of a SSO scheme.

The technique presented here uses the Ensemble Kalman Filter (EnKF) and Ensemble Kalman Smoother (EnKS) which are reviewed in detail in Evensen (2009). In this work, we do not use an augmented state to estimate parameters of GCMs, as is usually done for online estimation (Annan and Hargreaves, 2007; Ruiz et al., 2013); instead the state variables for the EnKF are only the physical parameters in this offline parameter estimation. As we do not have any knowledge of their temporal evolution, the state model is supposed to follow a random walk. In this way, we assume a non-negligible model error. An innovative part of our technique is that we also estimate the statistical parameters of the EnKF: (i) the covariance matrices of the Gaussian errors that control the weight of the state and the observation equations and (ii) the background state of the filter, typically an a priori knowledge of the physical parameters. Generally, these statistical parameters of the EnKF are prescribed values chosen by the user. In practice, this manual tuning does not ensure the filter convergence to the state of the system. To overcome this problem, the standard implementations of the EnKF use an inflation factor for the forecast and/or observational-error covariance matrices to avoid filter divergence. However, the main problem of this approach is the choice of the covariance inflation (additive or multiplicative) and the amplitude of the inflation. Several studies propose to estimate the inflation factors using the first moment estimation of the squared innovation (e.g. Wang and Bishop, 2003; Li et al., 2009; Liang et al., 2011), Bayesian approaches (e.g. Anderson, 2007; Miyoshi, 2011), or the second-order least-squares statistic of the squared innovation (Wu et al., 2013). The technique presented in this article does not need to use any inflation factor since the statistical parameters are non-deterministic values. Here, as the estimation is offline in a low-dimensional system, we estimate directly the entire error covariance matrices and the background state of the EnKF using a maximum likelihood approach. In particular, we use the iterative and efficient Expectation-Maximization (EM) algorithm introduced by Dempster et al. (1977). To our knowledge, the implemented technique in this work based on the combination of an ensemble Kalman filter with the EM algorithm has not been proposed previously in data assimilation.

The novel estimation technique is applied to the SSO scheme described in Lott and Miller (1997) and revised in Lott (1999). This SSO scheme computes the wind tendencies due to the subgrid-scale orography and is implemented in three GCMs: that of the Laboratoire de Météorologie Dynamique (LMDz), the ECHAM model which is the atmospheric component of the Earth System Model of the Max Planck Institute (MPI-ESM), and that of the European Centre for Medium-Range Weather Forecasts (ECMWF). It is known that weather forecast and climate models are sensitive to the physical parameters of SSO schemes (e.g. Palmer et al., 1986; Lott et al., 2005; Sigmond et al., 2008). Currently, this issue is still important since climate models now extend to the middle atmosphere where mountain gravity waves significantly affect the Brewer-Dobson circulation (McLandress and Shepherd, 2009). This circulation seems to intensify with climate change (Li et al., 2008). These results call for a re-evaluation of the SSO schemes in the middle-atmosphere-resolving models and in particular of the set of parameters used in the schemes. An optimization of the SSO schemes can help to evaluate better the potential effects of the orographic gravity wave drag on the westerlies in midlatitudes.

This article is organized as follows. First, we describe the SSO scheme and the datasets it uses in section 2. Then, in section 3, we present the statistical model used to estimate the physical parameters of the SSO scheme. The details of the estimation technique based on the EnKF, EnKS and the EM algorithm are explained in section 4. The estimation technique is applied to a column version (not a 3D version) of the subgrid-scale orography scheme. We then use two synthetic cases (i.e. without using real observations): an identical-twin experiment and a situation in which the horizontal resolution of the orography dataset is changed. We show the results in section 5. Conclusions are drawn and future work is outlined in section 6. In general, the unified notations of data assimilation given in Ide et al. (1997) are used here.

## 2. Data and model

### 2.1. General circulation model data

To conduct our offline estimation we used daily data from a simulation done with the LMDz GCM (Hourdin et al., 2006) using a horizontal resolution of 3.75° × 2.5° and 50 vertical levels with a model top at 5 hPa. We have extracted from this model the SSO scheme we want to optimize. To conduct the optimization, we limit ourselves to a one-month period, July 2000. The exact year itself is of little importance, since the run considered has a spin-up of several years, and was not constrained by forcings other than the sea-surface temperature and the land-sea ice cover. The particular month chosen is in midwinter in the Southern Hemisphere, when high wind speed conditions prevail over the southern Andes.

The SSO scheme we use represents mountain gravity wave drag and blocked flow drag following Lott and Miller (1997). It also introduces lateral lift to take into account the fact that narrow valleys are partially sheltered from the large-scale winds in the free troposphere (cf. Lott, 1999). The scheme was extended to the stratosphere in Lott et al. (2005) and this is the version we use. For completeness, the salient features of scheme are described here.

Before launching a simulation, subgrid-scale orography parameters are calculated in each model gridbox: the mountain minimum, mean, and maximum elevations, and the mountain departure from the mean is then characterized by its anisotropy, its orientation angle, its slope and its standard deviation. As we will see, when we change orography datasets these parameters change significantly and the most dramatic changes concern the evaluation of the slope. We will address these issues in section 5.2 and evaluate the changes to be done to the SSO scheme used in LMDz, when we make a transition from the 10 min of resolution US Navy orography dataset used in most current applications, to
<table><caption>Table 1. Physical parameters (defined in the text) of the SSO scheme, their assigned true values and their corresponding physical range.</caption><thead><tr><th>Physical parameter</th><th>True value</th><th>Range</th></tr></thead><tbody><tr><td>G</td><td>1</td><td>(0, 1.5)</td></tr><tr><td>C<sub>d</sub></td><td>1</td><td>(0, 1.5)</td></tr><tr><td>C<sub>l</sub></td><td>1</td><td>(0, 1.5)</td></tr><tr><td>H<sub>NC</sub></td><td>1</td><td>(0, 1.5)</td></tr><tr><td>β</td><td>0.5</td><td>(0, 1)</td></tr><tr><td>Ri<sub>c</sub></td><td>0.25</td><td>(0, 2)</td></tr></tbody></table>

a more refined 2 min of resolution dataset. At each time step, the SSO scheme uses the background flow conditions predicted by the model at a given gridpoint (i.e. the horizontal components of the winds, the temperature and the surface pressure), and predicts the effect of the SSO on the large-scale flow at all model levels.

The SSO scheme uses a set of six non-dimensional parameters of order O(1) which characterize the mesoscale and synoptic-scale effects of the mountain on the large-scale flow. The first three parameters directly scale the forces associated with the different processes parameterized: the gravity wave drag G, the low-level blocked flow drag C<sub>d</sub>, and the low-level lift C<sub>l</sub> that enhances large-scale vortex compression to represent valley sheltering. The other three parameters used are the low-level flow blocking depth H<sub>NC</sub>, the fraction β of the gravity wave drag that propagates toward the free troposphere and aloft, and the critical Richardson number Ri<sub>c</sub> that is used to predict when the mountain waves break. In Table 1, we give the values of each parameter used in operation, and also the range of values we will consider as plausible when we retune the SSO scheme.

The profiles of wind tendencies given by the SSO scheme are very sensitive to the value of the six non-dimensional parameters as in other schemes. The parameters used in the past were motivated by decades of research on mountain flow dynamics, and by a few experimental campaigns conducted over specific areas like the Pyrenees in France (e.g. Lott, 1995, gives a motivation for the lift based on the PYREX campaign). Although satellite data combined with high-resolution simulations could also be used in the future (Hertzog et al., 2012), it remains the case that local tunings will probably still be needed, at least near places where the drag forces can potentially be very important. For this reason, and also because the methodology we propose is well adapted to handle one-column models (Posselt and Bishop, 2012), this is the strategy we have followed in the present work, where we imagine that an observational campaign takes place near the Perito Moreno Glacier in the Andes (46°S, 71°W) (dot in Figure 1). There, the mountains are characterized by an important anisotropic shape and strong variation in altitude (the standard deviation is 295 m for a mean altitude of 531 m and a peak of 1513 m). These topographical conditions represented in Figure 1(a) are ideal to study mountain-induced forces especially in high surface wind speed conditions such as in Figure 1(b). Indeed, this geographical location gave one of the largest subgrid-scale mountain drag amplitudes on the Earth in a preliminary spatial analysis for July 2000 in which we computed globally the subgrid-scale mountain drag with the scheme.

## 2.2. Preliminary tests

To give a preliminary view of the scheme outputs, Figure 2 shows the tendencies predicted by the SSO scheme,

$$y(t_k) = \mathcal{F}\{\theta, \mathbf{Z}(t_k)\}, \quad (1)$$

where $\theta = (H_{NC}, C_d, Ri_c, G, C_l, \beta)$, and $\mathbf{Z}(t_k)$ is a generic notation for the vertical profiles of the horizontal winds and temperature. In Eq. (1), the vector $y(t_k)$ has $m = 100$ values each of the 31 days $k \in \{1, \dots, K = 31\}$ of July 2000. Each day, the first 50 values correspond to the zonal tendencies at the 50 model levels, and the last 50 values to the meridional tendencies.

From a preliminary temporal analysis during July 2000 at the chosen location, we distinguish two characteristic regimes of wind profiles in terms of the resulting induced SSO tendencies. Examples of these two regimes are shown in Figure 3. On 5 July 2000 (dotted line), the wind profile shows low and constant wind speeds with altitude. On 25 July (dashed line) the profile shows higher surface wind speeds and an increase of the wind with height, due to the presence of the subpolar jet in the region. The zonal and meridional components of the SSO tendency for the two cases are shown in Figure 2. The free physical parameters of the scheme are set to $H_{NC}^t = 1$, $C_d^t = 1$, $Ri_c^t = 0.25$, $G^t = 1$, $C_l^t = 1$ and $\beta^t = 0.5$. These set of ‘true’ parameter values $\theta^t$ were proposed by Lott (1999) and generate the ‘true’ tendency denoted $y^t(t_k)$. From Figure 2, on 25 July 2000, we find large tendencies whereas on 5 July 2000 the effect is much weaker due to the low wind speed conditions. The 95th percentile envelope around the mean value for the month of July 2000 is also shown. It indicates that the predicted SSO tendencies tend to be small in the mid-levels and larger at levels corresponding to the peak of the mountain (800 hPa) and to the tropospheric jet (250 hPa).

To evaluate how the outputs of the scheme vary with the different parameters, we evaluate the cost function

$$J(t_k) = \{y^t(t_k) - y(t_k)\}^\top \{y^t(t_k) - y(t_k)\}, \quad (2)$$

where the transpose notation $^\top$ is used, so that the square differences are summed over all altitudes and over the two components (zonal and meridional). Firstly, the cost function given in Eq. (2) is computed by changing one physical parameter and by fixing the other parameters to the true values.

Figure 4(a) shows the sensitivity of J as a function of H<sub>NC</sub> and C<sub>l</sub> for the state found on 5 July 2000 which is a situation with low surface winds as shown in Figure 3. The cost function associated with H<sub>NC</sub> shows a non-quadratic behaviour, representing a nonlinear sensitivity in the derivative of J. The parameter C<sub>l</sub> shows a quadratic cost function so that its sensitivity is linear. The other four parameters also show a linear sensitivity, as found for the C<sub>l</sub> parameter, so that their cost functions are not shown. Figure 4(b) shows the sensitivity of J in high surface wind speed conditions, on 25 July 2000. The sensitivity of J to H<sub>NC</sub> is enhanced by a factor of 10<sup>3</sup> in strong wind speed conditions compared to the sensitivity for the weak wind speed case. Finally, a relatively weaker enhancement of the J sensitivity is found for high surface wind speed conditions to certain physical parameters, e.g. C<sub>l</sub> (squares) and C<sub>d</sub> (not shown here) compared to the enhancement of H<sub>NC</sub> sensitivity between low and high surface wind conditions. In Figure 4(b), a saturation of the cost function is found close to the global minimum for H<sub>NC</sub> > 1.1 (circles) in high wind speed conditions. This behaviour can be explained as follows. As the surface wind increases and H<sub>NC</sub> increases, the blocked flow depth decreases and eventually reaches 0 (Eqs. (4) and (9) in Lott and Miller, 1997). At this point, the parameter H<sub>NC</sub> becomes saturated since an increase of its value cannot change the blocked flow depth to negative values. Therefore, H<sub>NC</sub> values larger than this critical value cannot affect the SSO predictions.

In a second sensitivity experiment, the cost function given in Eq. (2) is computed changing two physical parameters simultaneously. Figure 5 shows the cost function as a function of the parameters H<sub>NC</sub> and G. H<sub>NC</sub> is correlated with G. The ten smallest values of the cost function are indicated by black dots in Figure 5. They underline the fact that the global minimum region of the cost function (intersection of the two dashed black lines) is not well defined. On the contrary, in Figure 5(a), a large region of very low sensitivity close to the global minimum is highlighted, especially in low wind speed conditions where the sensitivity of J is reduced. In this region of the cost function, there is a negative correlation between H<sub>NC</sub> and G.
<figure><img src="image_58.png" /><figcaption>Figure 1. (a) Topography of the south Andes and (b) surface winds on 25 July 2000. The location of the chosen mountain peaks is close to the Perito Moreno Glacier (46°S, 71°W) and represented by a dot.</figcaption></figure>

### 3. Nonlinear Gaussian state-space model

To estimate the $n = 6$ physical parameters in $\theta$ via our filtering technique, we first need to make them stochastic. We denote them as $\mathbf{x}$ and we say it is the 'state of the system'. The state evolution is given by a Gaussian random walk,

$$ \mathbf{x}(t_k) = \mathbf{x}(t_{k-1}) + \eta(t_k), \quad (3) $$

where the $n$-dimensional stochastic random vector $\{\eta(t_k)\}_{k \in \{1,...,K\}}$ represents an additive perturbation at each time $t_k$. We assume that the perturbations are Gaussiany distributed with zero mean and a constant in time $n \times n$ covariance matrix $\mathbf{Q}$. Equation (3) is taken as the state equation in our state-space model.

If we use directly the physical parameters $\theta$ as the state of the system, they can easily become negative or reach very large values, whereas the parameters in the SSO scheme are assumed to be always positive and of the order of unity. For this reason, we map the physical parameters $\theta$ on $\mathbf{x}$ by using the Gauss error function $\theta = \mathcal{G}(\mathbf{x})$, as sometimes used in data assimilation (Hu et al., 2010).

At the initial time of Eq. (3), we introduce an a priori knowledge of the physical parameters. We assume that this background information follows a Gaussian distribution given by the $n$-dimensional vector mean $\mathbf{x}^b$ and the $n \times n$ covariance matrix $\mathbf{B}$.

At time $t_k$, the zonal and meridional SSO tendencies given in Eq. (1) are assumed to be observed. They are stored in the $m$-dimensional stochastic random vector $\{\mathbf{y}(t_k)\}_{k \in \{1,...,K\}}$. The state vector at time $t_k$ is related to the observation by means of the observation equation defined by

$$ \mathbf{y}(t_k) = \mathcal{H}_k \{\mathbf{x}(t_k)\} + \mathbf{\epsilon}(t_k), \quad (4) $$

where the observation operator $\mathcal{H}_k$ is the nonlinear function defined by

$$ \mathcal{H}_k \{\mathbf{x}(t_k)\} = \mathcal{F}[\mathcal{G}\{\mathbf{x}(t_k)\}, \mathbf{Z}(t_k)], \quad (5) $$

where $\mathcal{F}$ is the SSO scheme Eq. (1) and $\mathcal{G}$ is the Gauss error function. In Eq. (4), we suppose that the $m$-dimensional stochastic random vector $\{\mathbf{\epsilon}(t_k)\}_{k \in \{1,...,K\}}$ is an additive zero-mean Gaussian error. The $m \times m$ covariance matrix of $\mathbf{\epsilon}(t_k)$ is denoted by $\mathbf{R}(t_k)$.
<figure><img src="image_59.png" /><figcaption>**Figure 2.** Vertical profiles of (a) zonal and (b) meridional tendencies generated by the SSO scheme at location 46°S, 71°W within the Andes. The grey shading denotes the 95th percentile envelope over the month of July 2000. The dotted and dashed lines correspond to the mountain drag profiles on 5 and 25 July 2000 respectively.</figcaption></figure>

<figure><img src="image_60.png" /><figcaption>**Figure 3.** Wind speed profiles from the LMDz GCM model on 5 (dotted line) and 25 (dashed line) July 2000 at location 46°S, 71°W within the Andes.</figcaption></figure>

As the sensitivity of *J* varies with the atmospheric conditions Z(*t*<sub>k</sub>), particularly with the surface wind speed, **R**(*t*<sub>k</sub>) is assumed in principle to vary with time.

The statistical parameters correspond to the vector and matrices that define the system (3) and (4). They are denoted by $\psi$. We use the term 'statistical parameters' of the state-space statistical model to distinguish from the six 'physical parameters' $\theta$ of the SSO scheme. The statistical parameters are the a priori

<figure><img src="image_61.png" /><figcaption>**Figure 4.** Cost function (Eq. 2)) as a function of the physical parameters $H_{NC}$ (circles) and $C_1$ (squares). The true values of the physical parameters are $H_{NC}^t = 1$ and $C_1^t = 1$. The results are given for the location 46°S, 71°W on (a) 5 July 2000 representing low surface wind speed conditions and (b) 25 July 2000 representing high surface wind speed conditions.</figcaption></figure>

probability density functions (PDFs) of the physical parameters, given by $\mathbf{x}^b$ and $\mathbf{B}$, and the covariance error matrices $\mathbf{Q}$ and $\mathbf{R}(t_k) \forall k \in \{1, ..., K\}$. We write $\psi = (\mathbf{x}^b, \mathbf{B}, \mathbf{Q}, \mathbf{R})$. The statistical parameters define the uncertainty of the state-space statistical model and play a central role on the quality and rate of convergence in the estimation of the physical parameters with the filtering and smoothing techniques described below.

The estimation of the statistical parameters $\psi$ is conducted maximizing the total likelihood function $\mathcal{L}$. This function is based on the PDF of the initial state $p\{\mathbf{x}(t_1)\}$, the conditional state evolution $p\{\mathbf{x}(t_k)|\mathbf{x}(t_{k-1})\}$ and the observations conditionally to the state $p\{\mathbf{y}(t_k)|\mathbf{x}(t_k)\}$. The three PDFs are assumed to be normally distributed with the respective mean and covariances: $\mathbf{x}(t_1) - \mathbf{x}^b$ and $\mathbf{B}$, $\mathbf{x}(t_k) - \mathbf{x}(t_{k-1})$ and $\mathbf{Q}$, $\mathbf{y}(t_k) - \mathcal{H}_k\{\mathbf{x}(t_k)\}$ and $\mathbf{R}(t_k)$. Finally, using the Markov property of the state-space model, the total likelihood function is the product of the PDF for all times $K$. It is given by

$$ \mathcal{L}(\mathbf{x}, \psi) = p\{\mathbf{x}(t_1)\} \prod_{k=2}^{K} p\{\mathbf{x}(t_k)|\mathbf{x}(t_{k-1})\} \times \prod_{k=1}^{K} p\{\mathbf{y}(t_k)|\mathbf{x}(t_k)\}. \quad (6) $$

In practice, this total likelihood function is approximated by its expectation conditionally to all the observations
<figure><img src="image_62.png" /><figcaption>Figure 5. Cost function (Eq. (2)) as a function of the physical parameters $H_{NC}$ (x-axis) and G (y-axis). The intersection of dashed black lines corresponds to the true physical parameters $H_{NC}^1 = 1$ and $G^1 = 1$. The black dots denote the physical parameters which generate the ten lowest values of the cost function. The results are given for the location 46°S, 71°W on (a) 5 July 2000 representing low surface wind speed conditions and (b) 25 July 2000 representing high surface wind speed conditions.</figcaption></figure>

$y_{1:K} = y(t_1), \dots, y(t_K)$. This requires the computation of the state smoothed probabilities to be described below.

## 4. Estimation technique

The algorithm to estimate the physical and statistical parameters is described concisely here. A diagram with the main steps of the algorithm is shown in Figure 6. The algorithm starts with a proposed set of statistical parameters $\hat{\psi}^{(1)}$ which do not need to be known precisely. Then, the statistical parameters are estimated by maximizing the total likelihood function using the EM algorithm. A loop is initiated which is composed by an expectation and a maximization step. The expectation step computes the expectations given in Appendix A via the EnKS. The maximization step consists basically in computing the optimal $\hat{\psi}^{(j)}$ from the known analytical expressions given in Appendix B. At each iteration j of the EM algorithm, we compute the innovation likelihood l given in Appendix C. It is commonly used to evaluate the quality of the state estimates and to compare state-space models with different statistical parameters (Cappé et al., 2005, p. 140 give more details). If the innovation likelihood does not change significantly, the last estimated $\hat{\psi}^{(j)}$ is returned. These optimal statistical parameters given by the EM algorithm are finally used to initiate a last EnKF run which estimates the physical parameters.

<figure><img src="image_63.png" /><figcaption>Figure 6. Diagram of the method based on the maximum likelihood estimates of the state-space model Eqs (3) and (4).</figcaption></figure>

### 4.1. Expectation-maximization algorithm

The maximum likelihood estimates of the statistical parameters $\psi$ are conducted using the EM algorithm proposed by Dempster et al. (1977). This is a classical method used in the case of incomplete or missing data. This iterative algorithm is based on two steps: the expectation of the total log-likelihood function (E step) and its maximization with respect to $\psi$ (M step). The EM algorithm begins with an initial set of statistical parameters $\hat{\psi}^{(1)}$. Then, repeating the E and M steps, the sequence of estimates $\hat{\psi}^{(j)}$ yields increasing values of the expected log-likelihood and converges to the maximum likelihood estimates.

At iteration j, the E step consists of computing the expected total log-likelihood function conditionally to the total observations and the previously estimated statistical parameters. It is given by

$$Q(\psi | \hat{\psi}^{(j-1)}) = E \left[ \log \{ \mathcal{L}(\mathbf{x}, \psi) \} | y_{1:K}, \hat{\psi}^{(j-1)} \right]. \quad (7)$$
In the case of nonlinear state-space statistical models, the exact smoothed probabilities are not computable. Thus, we use the Monte Carlo approximations given by the EnKS. The conditional expectations given in Appendix A are then computed.

The M step consists of maximizing $Q(\psi | \hat{\psi}^{(j-1)})$ with respect to $\psi$. We obtain a direct analytic form of the maximum likelihood estimates. The expressions are given in Appendix B. The derivations are not presented here (cf. Tandeo et al., 2011, for more details).

## 4.2. Ensemble Kalman filter

The EnKF algorithm used here is an adaptation of the one proposed by Burgers et al. (1998).

In the initial step of the EnKF algorithm, at time $t_1$, an ensemble of $\mathbf{x}_i$ composed of $N$ members is randomly generated. The members of the ensemble follow a Gaussian distribution given by the vector mean $\mathbf{x}^b$ and the covariance matrix $\mathbf{B}$. The $N$ initial members are stored in the vectors $\mathbf{x}_i^f(t_1) \forall i \in \{1, ..., N\}$.

In the update step, at each time $t_k$, we randomly generate $N$ samples of $\eta_i$ and $\epsilon_i \forall i \in \{1, ..., N\}$ with respective covariances $\mathbf{Q}$ and $\mathbf{R}(t_k)$. Then, following Eq. (3), the $i$-member of the updated state is given by

$$ \mathbf{x}_i^f(t_k) = \mathbf{x}_i^a(t_{k-1}) + \eta_i(t_k), \quad (8) $$

and the mapping from the forecast state space to the observational space of the $i$-member is computed as

$$ \mathbf{y}_i^f(t_k) = \mathcal{H}_k \{ \mathbf{x}_i^f(t_k) \}. \quad (9) $$

The $N$ members of the ensemble are used to estimate the sample means of the propagated state in the state space and in the observational space denoted by $\mathbf{x}^f(t_k)$ and $\mathbf{y}^f(t_k)$ respectively.

In the analysis step, we follow the Pham (2001) methodology which avoids the linearization of the observational operator. The Kalman gain is computed with

$$ \mathbf{K}(t_k) = \mathbf{P}_{xy}^f(t_k) \{ \mathbf{P}_{yy}^f(t_k) + \mathbf{R}(t_k) \}^{-1}, \quad (10) $$

where $\mathbf{P}_{xy}^f(t_k)$ is the sample cross-covariance matrix and $\mathbf{P}_{yy}^f(t_k)$ is the sample covariance matrix, which are determined by

$$ \mathbf{P}_{xy}^f(t_k) = \frac{1}{N-1} \sum_{i=1}^{N} \{\mathbf{x}_i^f(t_k) - \mathbf{x}^f(t_k)\} \{\mathbf{y}_i^f(t_k) - \mathbf{y}^f(t_k)\}^{\top} \quad (11) $$

and

$$ \mathbf{P}_{yy}^f(t_k) = \frac{1}{N-1} \sum_{i=1}^{N} \{\mathbf{y}_i^f(t_k) - \mathbf{y}^f(t_k)\} \{\mathbf{y}_i^f(t_k) - \mathbf{y}^f(t_k)\}^{\top}. \quad (12) $$

In our case, the number of observations is larger than the dimension of the state space ($m > n$), so that the matrix $\mathbf{P}_{yy}^f(t_k) + \mathbf{R}(t_k)$ is ill-conditioned making the matrix inversion difficult. Therefore, as described in Evensen (2009), Chapter 14, we compute the pseudo-inverse of $\mathbf{P}_{yy}^f(t_k) + \mathbf{R}(t_k)$, taking into account 99% of the information given by the eigenvalues. Having $\mathbf{K}(t_k)$ from Eq. (10), the $N$ members of the ensemble are then updated by

$$ \mathbf{x}_i^a(t_k) = \mathbf{x}_i^f(t_k) + \mathbf{K}(t_k)\mathbf{d}_i(t_k), \quad (13) $$

where the $m$-dimensional $\mathbf{d}_i(t_k) \forall i \in \{1, ..., N\}$ are the $N$ innovation vectors in which we use perturbed observations such as $\mathbf{d}_i(t_k) = \mathbf{y}(t_k) + \epsilon_i(t_k) - \mathbf{y}_i^f(t_k)$. Note that the sample covariance of the $N$ innovations is $\mathbf{P}_{yy}^f(t_k) + \mathbf{R}(t_k)$. Finally, the updated analyzed state is represented by the sample mean $\mathbf{x}^a(t_k)$ and the sample covariance $\mathbf{P}^a(t_k)$.

## 4.3. Ensemble Kalman smoother

The backward recursions correspond to the EnKS algorithm proposed by Evensen and Van Leeuwen (2000). It uses the results of the EnKF computed above.

In the initial step of the EnKS algorithm, at time $t_K$, we use the members of the filtered state, $\forall i \in \{1, ..., N\}$, such as $\mathbf{x}_i^s(t_K) = \mathbf{x}_i^a(t_K)$ and $\mathbf{P}^s(t_K) = \mathbf{P}^a(t_K)$.

Then, we proceed backward from $k = K - 1$ to $k = 1$. At each time $t_k$, we compute

$$ \mathbf{x}_i^s(t_k) = \mathbf{x}_i^a(t_k) + \mathbf{K}^s(t_k) \{ \mathbf{x}_i^s(t_{k+1}) - \mathbf{x}_i^f(t_{k+1}) \}, \quad (14) $$

where $\mathbf{K}^s(t_k)$ is the $n \times n$ Kalman smoother gain matrix given by $\mathbf{P}^a(t_k) \{ \mathbf{P}^f(t_{k+1}) \}^{-1}$. The Gaussian distribution of the updated state estimate is given by the sample mean and covariance respectively denoted by $\mathbf{x}^s(t_k)$ and $\mathbf{P}^s(t_k)$. The sample covariance of the state between two consecutive times is computed using

$$ \mathbf{P}^s(t_k, t_{k-1}) = \frac{1}{N-1} \sum_{i=1}^{N} \{\mathbf{x}_i^s(t_k) - \mathbf{x}^s(t_k)\} \{\mathbf{x}_i^s(t_{k-1}) - \mathbf{x}^s(t_{k-1})\}^{\top}. \quad (15) $$

# 5. Results

## 5.1. Identical-twin experiment

In order to evaluate the technique, twin experiments are used. In this case, the observations are obtained under the assumption of a perfect model; in other words, the SSO scheme is assumed to give the true tendencies when the physical parameters $\theta^t = (1, 1, 0.25, 1, 1, 0.5)$ are used as the true parameters. Then we suppose that the state, i.e. the physical parameters, is unknown and we try to estimate it via the state-space model (3) and (4) using the generated mountain drag observations. As schematized in Figure 6, we estimate the statistical parameters $\psi = (\mathbf{x}^b, \mathbf{B}, \mathbf{Q}, \mathbf{R})$ of the state-space model via the EM algorithm in order to improve the estimation of the physical parameters, $\theta$. At iteration $j = 1$, we deliberately initialize the state vector $\hat{\mathbf{x}}_b^{(1)}$ far from the true state values (corresponding to the true parameters). The corresponding covariance $\hat{\mathbf{B}}^{(1)}$ is chosen as the unit matrix $\mathbf{I}_6$ to generate large initial spreads of the members. Throughout the filter evolution, the members are randomly perturbed by the constant covariance matrix $\hat{\mathbf{Q}}^{(1)} = 0.1 \times \mathbf{I}_6$ in Eq. (8). The covariance of the measurement errors in Eq. (9) is set to $\hat{\mathbf{R}}^{(1)}(t_k) = 1000 \times \mathbf{I}_{100} \forall k \in \{1, ..., K\}$, which is of the same order as the mean value of the cost function $J$ given in Eq. (2). We use $N = 100$ members and 25 iterations of the EM algorithm.

The innovation log-likelihood function and the total RMSE of the physical parameters for the conducted twin experiments are shown in Figure 7(a) as a function of the EM iteration. The results indicate that the innovation log-likelihood is a good synthetic indicator of the filter quality which follows the inverse variations of the total RMSE. In Figure 7(b), we decompose the total RMSE for each physical parameter. We find a good convergence of all the physical parameters after $j = 10$ iterations except for $C_d$, $H_{NC}$ and $G$, which need more EM iterations. The evolution of these two last physical parameters as a function of time for different iterations ($j = 1, 10, 25$) of the EM algorithm are shown in Figure 8. For both physical parameters, the EM algorithm is able to adapt the filter conditions and to give, along the iterations $j$, more and more accurate initial distributions of the physical parameters (given by the $\mathbf{x}^b$ and $\mathbf{B}$ maximum likelihood estimates). However, at the last iteration $j = 25$, the temporal convergence (near $k = 20$) is higher than the other physical parameters (not shown here). Note that the results using deterministic values of $\psi$, instead of estimating them via the maximum likelihood method, show the inability of the filter to converge to the solution $\theta^t$. This is shown with the
<figure><img src="image_64.png" /></figure>

<figure><img src="image_65.png" /><figcaption>**Figure 7.** Evolution of (a) the innovation log-likelihood (left y-axis; blue line), the total physical parameter RMSE summed by time (right y-axis; red line) and (b) detailed RMSE of each physical parameter over *j* = 25 iterations of the EM algorithm.</figcaption></figure>

blue curves of Figure 8 corresponding to the first iteration of the EM algorithm (i.e. this could be interpreted as a standard EnKF estimation). Even if we use more realistic but uniform values of covariance matrices $\hat{\mathbf{B}}^{(1)}$, $\hat{\mathbf{Q}}^{(1)}$ and $\hat{\mathbf{R}}^{(1)}(t_k) \forall k \in \{1, \dots, K\}$, the standard EnKF is unable to converge to a stable and accurate solution.

Figure 9(a) shows the matrix **Q** after 25 EM iterations. A negative correlation between the $H_{NC}$ and $G$ physical parameters is clearly detected. This confirms the observation we made from Figure 5 in the weak sensitivity region of the cost function $J$. The elements of the **Q** estimated by the maximum likelihood method for the covariance between $H_{NC}$ and $G$ and the variances of $H_{NC}$ and $G$ are respectively $-1.5 \times 10^{-5}$, $1.5 \times 10^{-5}$ and $3 \times 10^{-5}$. These variances correspond to the optimal perturbations of the members in Eq. (3) at each time of the filter. Note that the amplitude of **Q** tends to decrease with the iterations of the EM algorithm since the model becomes perfect and the observations are produced with the optimal physical parameters. Concerning the estimated amplitude of the observation-error covariance $\mathbf{R}(t_k)$, i.e. the covariance of $\epsilon(t_k) \forall k \in \{1, \dots, K\}$, it varies with the forcing terms, particularly the surface wind speed conditions. The results for the low and high wind speed conditions are shown in Figure 9(b) and (c). We distinguish different parts on these estimated matrices. The top left and the bottom right parts correspond respectively to the zonal and meridional error covariances of the observation equation given in Eq. (4). The top right and bottom left parts correspond to the cross-covariance between the zonal and meridional components. The x- and y-axes indicate the vertical level of the different components. For instance, the level 1000 hPa is given by the indices 1 and 51 whereas the level 5 hPa is given by the indices 50 and 100. The results indicate a checkerboard structure in the covariances inside groups of vertical levels and especially a larger variability

<figure><img src="image_66.png" /></figure>

<figure><img src="image_67.png" /><figcaption>**Figure 8.** Evolution of the (a) $H_{NC}$ and (b) $G$ estimates with their 95% confidence intervals over time for different iterations of the EM algorithm: $j = 1$ (blue), $j = 10$ (red) and $j = 25$ (black). The straight lines correspond to the true physical parameter values. The physical parameters ($\theta$, not **x**) are shown.</figcaption></figure>

of the observation error in the levels close to the surface in both cases. We also observe a cross-covariance between the zonal and meridional error terms at this altitude for the low wind speed conditions. The main difference between the two estimated matrices is the amplitude of the variability. In strong wind speed conditions, the variance is globally enhanced by a factor of $10^3$. Therefore, the method proposed here is able to model a flow-dependent (typically the wind speed) and not necessarily diagonal error covariance matrix $\mathbf{R}$. Miyoshi et al. (2012) have also proposed to retrieve the shape of $\mathbf{R}$ in a data assimilation problem conducting twin experiments. More precisely, they extended the adaptive estimation method proposed by Li et al. (2009) to include off-diagonal terms of $\mathbf{R}$.

We make two comments on results that are not shown here. Firstly, the use of $N = 500, 1000$ members in the ensemble (not shown) gives similar results as the case with $N = 100$ presented here. Thus, an ensemble of 100 members is sufficient to capture the highly nonlinear behaviour of the SSO scheme and to estimate properly the statistical parameters of the state-space system. Secondly, the maximum likelihood statistical parameter $\psi$ estimates are independent of the initial conditions of the EM algorithm. Different initial guess parameters $\hat{\mathbf{x}}_b^{(1)}$ and different covariances $\hat{\mathbf{B}}^{(1)}$ give similar rates of convergence.

## 5.2. *Changes in orography resolution*

When the resolution of a GCM is increased, or when a new dataset is used to feed the physical parametrizations, the physical parameters of the GCM need to be adjusted. There is no systematic
<figure><img src="image_68.png" /><figcaption>**Figure 9.** Maximum likelihood estimates after *j* = 25 iterations of the EM algorithm of (a) **Q**, (b) **R**(*t*<sub>5</sub>) and (c) **R**(*t*<sub>25</sub>). The straight lines in (b) and (c) denote the limit between the zonal and meridional mountain drag error covariance of the observation equation.</figcaption></figure>

way to produce these adjustments in the schemes so far. The
technique introduced in this work can be used to do this. In
particular, the standard parameters that are currently used in the
SSO scheme shown in Table 1 have been manually tuned using
PYREX data by Lott and Miller (1997). This set of parameters are
used operationally in the LMDz model. The tuning was conducted
with a version of the SSO scheme that uses the low-resolution
orography (10' × 10'; Figure 10(a)). Suppose that the higher
resolution (2' × 2'; Figure 10(b) and NOAA, 2001) orography
dataset is used to improve LMDz at a given horizontal resolution.
The parameters of the scheme should be adjusted for this new
orography dataset. We conducted an experiment to examine if

the technique is able to determine a new set of optimal parameters
for this high-resolution orography dataset. The conducted data
assimilation experiments use the SSO tendencies predicted by the
SSO scheme using the low-resolution orography as observations.
The assimilation then uses the SSO tendencies predicted with a
higher-resolution orography dataset. In this way, the operator *H*<sub>*k*</sub>
used in the assimilation has an error.

In this experiment in which the SSO scheme is not ‘perfect’
due to the resolution change, we found that the results depend on
the initial guess conditions, in opposition to the identical-twin
experiments. As convergence could not be reached easily, one
hundred filter experiments with different random initial guess
conditions $\hat{\mathbf{x}}_b^{(1)}$ of the EM algorithm were conducted. Among
these 100 experiments, we find results with recurrent estimations
reach the same log-likelihood as those shown in Figure 11(a). The
parameter estimations after $j = 25$ EM iterations for five selected
cases are shown in Figure 11(b).

From Figure 11(b), we notice that there is one parameter
which does not need to be changed much when the resolution
is changed –$C_j$. This is not a surprise since $C_j$ is an almost
linear lift coefficient, which is related to a mountain lift force
whose amplitude varies linearly with the difference between the
mountain and valley height. We also find that $C_d$ needs to
be reduced by a factor ≈ 2. Considering Eq. (16) in Lott and
Miller (1997), the scheme measures the number of mountains
in a subgrid-scale area, and multiplies the low-level drag by
this number of ridges. This yields a multiplicative factor in the
mountain slope. When we move to a higher-resolution grid, the
estimate of the slope necessarily increases, so $C_d$ needs to decrease.
The same conclusion could be drawn for the parameter *G* which
controls the gravity wave drag, but here the technique gives two
possible solutions. One where *G* is almost unchanged or has
a weak increase, and one where it is decreased substantially, as
expected. As the solution with unchanged *G* is the most surprising,
it is important to notice that this is also related to a smaller *β*;
they therefore correspond to more trapped waves which apply
more low-level drag. As at the low level, it is $C_d$ that essentially
controls the drag. We have therefore increased the gravity wave
drag by increasing *G* but placed that drag at low level where the
effect is small compared to that of $C_d$. Another important result
of the analysis is that the value of the critical Richardson number
clearly converges to $Ri_c = 1.5$. As this high-resolution orography
case likely has larger-amplitude gravity waves, this larger
Richardson number than the one used with the low-resolution
orography dataset needs to be enhanced so that the waves
propagate at high levels without breaking systematically at lower
levels.

In general, the parameter estimations, except for *Ri*<sub>c</sub>, present
a very large spread, particularly for those parameters acting at
low levels. For these, it should be remembered that the drag
at low levels is always treated via implicit methods in part for
stability, and in part because overestimated drags could yield
wind reversals at low levels, which contradict the nature of drag
forces. Clearly, the assimilation technique indicates that some
physical considerations should be given to make these parameters
more efficient in controlling the drag. Among the possibilities,
the SSO scheme does not consider that, when there are several
mountains in a gridbox area, some sheltering should be taken into
account not to decelerate the same flow twice in succession. This
is currently handled implicitly by the scheme, but the low-level
drag should take into account this horizontal sheltering when we
increase the orography resolution.

Figure 11(b) shows that the filter converges towards two
possible optimal states, in which $H_{NC}$ and *G* clearly present
bimodal distributions. This result is associated with the high
correlation that was found in the cost function between $H_{NC}$ and
*G* (Figure 5). The presence of model error in this imperfect model
experiment appears to add complexity to the cost function with
the presence of these two local minima. This is consistent with the
results obtained by Schirber *et al.* (2013) in an online parameter
<figure><img src="image_69.png" /><figcaption>**Figure 10.** (a) Low 10′ × 10′ and (b) high 2′ × 2′ topographical resolution near location 46°S, 71°W (black dot) in the southern Andes.</figcaption></figure>

estimation under the presence of model error. From a physical perspective, this bimodal result is not surprising since large $H_{NC}$ yields low blocking levels, and a more efficient mountain elevation to excite gravity waves. When there is a larger amount of gravity wave drag, a good fraction of the corresponding gravity waves is likely to break at low level; this may be an effect hidden in the low-level drag discussed earlier. Also, this bimodality may be inherent in the nonlinear low-level flow dynamics the scheme tries to represent.

Figure 12 shows the five profiles of the SSO tendency intensity (i.e. the norm of the SSO tendency) generated with the estimated parameters for weak (5 July 2000) and strong (25 July 2000) surface wind conditions. In both surface wind conditions, the sets of estimated parameters with large *G* tend to underestimate the low-level drag (between 900 and 1000 hPa) and to overestimate the drag at high levels (between 650 and 900 hPa). On the other hand, the sets of estimated parameters with small *G* (and large $H_{NC}$) tend to overestimate the drag at low levels and also at higher levels (but they are relatively closer to the observed than in the cases with large *G* at those levels). The spread in the *β* parameter also appears to play a role.

## 6. Conclusion and outlook

In this article, we use a filtering technique to estimate the physical parameters of a subgrid-scale orographic scheme. The estimation is conducted offline, without estimating the state of the atmosphere and thus reducing the size of the state vector. As forcing terms, we use simulations of a general circulation model. The estimation problem is written as a nonlinear state-space system. This formulation is flexible and overcomes the main difficulties such as the boundaries on the physical parameters (strictly positive), the unknown background covariances and the high nonlinearity of the orographic scheme. In this state-space model, we suppose that the state and observation equations have additive Gaussian noise and that we know the *a priori* distribution of the physical parameters. The choice of these statistical parameters constitutes an important condition of convergence of the system to the true physical parameters. Thus, we estimate them via a maximum likelihood method. We use an iterative algorithm that computes the expected total log-likelihood function and maximizes it with respect to the statistical parameters.
<figure><img src="image_70.png" /><figcaption>**Figure 11.** Evolution of five cases, with $N = 100$ members and the high-resolution orographic scheme, of (a) the innovation log-likelihood along the EM iterations and (b) the physical parameters estimated by the EnKF over time at iteration $j = 25$ of the EM algorithm. The straight dashed lines in (b) denote the physical parameter values of the low-resolution orographic scheme. In (b), the physical parameters ($\theta$, not $x$) are shown.</figcaption></figure>

The estimation technique is evaluated in a single vertical column and using synthetic observations (i.e. without using real observations but those produced by the SSO scheme). We imagine that an observational campaign takes place near the Perito Moreno Glacier in the Andes, where the topographical conditions are ideal to study mountain drag. First, we use twin experiments: we prescribe a true set of physical parameters and generate synthetic observations of mountain drag. Then, we apply the estimation technique using these generated observations and compare the estimated parameters to the true ones. The results indicate a convergence of the filter to the true parameters after ~20 iterations of the EM algorithm. Even if the user initializes the error covariances and initial guess conditions with inappropriate values, these statistical parameters are iteratively updated and will converge towards the optimal values. The technique is able to detect correlations between parameters, to weight the observations as a function of the external forcing terms and to generate adaptive *a priori* information on the parameters. This overcomes the results obtained with deterministic values of statistical parameters which are usually arbitrarily prescribed since they are unknown.

We also examined whether the estimation technique is useful to determine whether the physical parameter should be changed when the horizontal resolution of an input dataset of the general circulation model is increased. In this case, the SSO scheme is imperfect and our filter takes into account this model error adding Gaussian noises controlled by time-dependent covariance matrices. The results show that our technique is a useful tool to determine the changes in the parameter when the resolution

<figure><img src="image_71.png" /><figcaption>**Figure 12.** Profiles of the SSO tendency amplitude generated with the five sets of estimated parameters obtained with the SSO scheme using the high-resolution dataset on (a) 5 July 2000 and (b) 25 July 2000. The profiles used as observations and generated with the SSO scheme using the low-resolution dataset are shown by dashed lines.</figcaption></figure>

of the input orography dataset increases. However, model error degrades the estimated drag profiles; some features of the observed drag profile in Figure 12 cannot be reproduced by the estimated drag profiles that use the high-resolution orography dataset. A technique with model bias treatment as in Dee and Da Silva (1998) may be required to diminish the differences in the drag profiles. We also detected that some parameters may have a range of values for which the RMSE and the likelihood (cf. Figure 11(a)) almost do not change. These results show that there is no sensitivity to these parameters and therefore show that a precise value for these parameters is not important. We attribute this to the fact that, in the SSO scheme, a lot of drag is applied at low level and handled implicitly. In the scheme also, the low-level drag is multiplied by the number of ridges present in the gridbox area, a number which is around 1 or 2 for the US Navy $10' \times 10'$ dataset, and which becomes much larger when a more refined dataset is used. Ideally, we should take into account that, when a mountain exerts a drag, a wake downstream is associated with it, so that for mountains in the lee but still in the gridbox, the incident flow should be much reduced. Currently, when we increase the orography resolution, this effect is handled by an implicit treatment. Numerically this situation is satisfying, but clearly call for further understanding of the dynamical sheltering, and its impact on the large-scale flow. It may explain the difficulty in estimating the parameters and the difference between the drags generated with low- and high-resolution orographic datasets in Figure 12.
The technique presented here is an efficient method to resolve offline physical parameter estimation. The advantages are (i) the flexibility of the state-space formulation which can be applied to a large number of applications, (ii) the ability to estimate the background state and the eventually flow-dependent (and not necessarily diagonal) error covariance matrices Q and R of the EnKF, and (iii) the relative low computational cost of the technique where a relatively small number of members and few iterations of the EM algorithm are needed. One possible extension of the technique is the estimation of biases in more realistic cases with different kinds of model error.

This work is focused on the evaluation of the technique using first twin experiments and also an experiment with synthetic observations but using a higher-resolution orographic dataset so that the model (used in the data assimilation system) in this case is imperfect. In a real application, the technique requires vertical profiles of small-scale momentum forcing. We envisage two possible sources of this forcing which can be used to constrain orographic parameters. The most significant one is from intensive observational campaigns over mountainous areas. One of the most representative ones was PYREX (Bougeault et al., 1990). Currently, there are several proposed campaigns (over the Andes, over New Zealand and over Scandinavia) for intense measurements over mountains with aircraft, lidars and radiosondes. These combined instruments can give significant information on momentum fluxes and their divergences. These potential campaigns could be an important source of observational data to estimate parameters of the subgrid-orography schemes using the proposed technique. A second possible data source of small-scale momentum forcing can be obtained from data assimilation techniques. Pulido and Thuburn (2005) show that four-dimensional variational assimilation can be used to estimate the missing momentum forcing term in the model equations. The technique is applied to obtain missing momentum forcing profiles in the middle atmosphere; here a significant part of systematic model error can be associated with gravity wave drag since the other physical parametrization active at those levels—the radiative transfer scheme—contains well-known parameters. On the other hand, in the troposphere several parametrizations are coupled so that the source of missing momentum is not readily identifiable with a particular parametrization. Therefore, the data assimilation techniques might be potentially useful to constrain subgrid-orography schemes using only the momentum forcing profile in the stratosphere. However, the impact of model errors from different sources in the parameter estimation problem needs to be further investigated. Another point that needs to be further investigated in an actual application of this offline technique is the possible feedbacks between the parametrization and the low-level flow; these feedback processes can affect the optimal parameters.

Follow-up work could apply this technique for online parameter estimation in strongly nonlinear systems. A first step will be to evaluate the method in a low-dimension system. Parameter estimation in a low-dimensional model was previously done by Annan and Hargreaves (2004) using deterministic values of the background state and the error covariance matrices. The advantage of applying our technique is to estimate them properly via the EM algorithm. Some first simulations we have performed give promising results. A simplified version of the method may also be useful in a larger-dimension online parameter estimation problem, for instance when there are a few unknown statistical parameters which need to be estimated precisely.

## Acknowledgements

This work was supported by ANPCyT Argentina under grant PICT 2007 No. 411, the ECOS-Sud project DIAGAC and by the LEFE project IAC. The authors would like to thank Dr P. Ailliot and Dr J.J. Ruiz for their helpful comments.

# Appendix A

## E step

At each iteration j of the EM algorithm, we need the following conditional expectations:

$$ 
\begin{aligned}
E\left[\mathbf{x}(t_k)\mathbf{x}(t_{k-1})^{\top}|\mathbf{y}_{1:K}, \hat{\mathbf{\psi}}^{(j-1)}\right] \\
&= \mathbf{x}^{\mathrm{s}}(t_k)\mathbf{x}^{\mathrm{s}}(t_{k-1})^{\top} + \mathbf{P}^{\mathrm{s}}(t_k) \\
E\left[\mathbf{x}(t_k)\mathbf{x}(t_{k-1})^{\top}|\mathbf{y}_{1:K}, \hat{\mathbf{\psi}}^{(j-1)}\right] \\
&= \mathbf{x}^{\mathrm{s}}(t_k)\mathbf{x}^{\mathrm{s}}(t_{k-1})^{\top} + \mathbf{P}^{\mathrm{s}}(t_k, t_{k-1}) \\
E\left[\mathcal{H}_k\{\mathbf{x}(t_k)\}|\mathbf{y}_{1:K}, \hat{\mathbf{\psi}}^{(j-1)}\right] \\
&= \mathcal{H}_k\{\mathbf{x}^{\mathrm{s}}(t_k)\} \\
E\left[\mathcal{H}_k\{\mathbf{x}(t_k)\}\mathcal{H}_k\{\mathbf{x}(t_{k-1})\}^{\top}|\mathbf{y}_{1:K}, \hat{\mathbf{\psi}}^{(j-1)}\right] \\
&= \mathcal{H}_k\{\mathbf{x}^{\mathrm{s}}(t_k)\}\mathcal{H}_k\{\mathbf{x}^{\mathrm{s}}(t_{k-1})\}^{\top} + \mathbf{P}_{yy}^{\mathrm{s}}(t_k),
\end{aligned}
$$

where $\mathbf{P}_{yy}^{\mathrm{s}}(t_k)$ is the sample covariance of the $\mathcal{H}_k\{\mathbf{x}_i^{\mathrm{s}}(t_k)\} \forall i \in \{1, ..., N\}$.

# Appendix B

## M step

The maximum likelihood estimates of the statistical parameters are given by:

$$ 
\begin{aligned}
\hat{\mathbf{x}}_{\mathrm{b}}^{(j)} &= E\left[\mathbf{x}(t_1)|\mathbf{y}_{1:K}, \hat{\mathbf{\psi}}^{(j-1)}\right], \\
\hat{\mathbf{B}}^{(j)} &= \mathrm{Var}\left[\mathbf{x}(t_1)|\mathbf{y}_{1:K}, \hat{\mathbf{\psi}}^{(j-1)}\right], \\
\hat{\mathbf{Q}}^{(j)} &= \frac{1}{T-1}\sum_{k=2}^{T} E\left[\mathbf{x}(t_k)\mathbf{x}(t_{k-1})^{\top}|\mathbf{y}_{1:K}, \hat{\mathbf{\psi}}^{(j-1)}\right] \\
&\quad - \frac{1}{T-1}\sum_{k=2}^{T} E\left[\mathbf{x}(t_k)\mathbf{x}(t_{k-1})^{\top}|\mathbf{y}_{1:K}, \hat{\mathbf{\psi}}^{(j-1)}\right] \\
&\quad - \frac{1}{T-1}\sum_{k=2}^{T} E\left[\mathbf{x}(t_k)\mathbf{x}(t_{k-1})^{\top}|\mathbf{y}_{1:K}, \hat{\mathbf{\psi}}^{(j-1)}\right]^{\top} \\
&\quad + \frac{1}{T-1}\sum_{k=2}^{T} E\left[\mathbf{x}(t_{k-1})\mathbf{x}(t_{k-1})^{\top}|\mathbf{y}_{1:K}, \hat{\mathbf{\psi}}^{(j-1)}\right], \\
\hat{\mathbf{R}}^{(j)}(t_k) &= \mathbf{y}(t_k)\mathbf{y}(t_k)^{\top} \\
&\quad - E\left[\mathcal{H}_k\{\mathbf{x}(t_k)\}|\mathbf{y}_{1:K}, \hat{\mathbf{\psi}}^{(j-1)}\right]\mathbf{y}(t_k)^{\top} \\
&\quad - \mathbf{y}(t_k)E\left[\mathcal{H}_k\{\mathbf{x}(t_k)\}|\mathbf{y}_{1:K}, \hat{\mathbf{\psi}}^{(j-1)}\right]^{\top} \\
&\quad + E\left[\mathcal{H}_k\{\mathbf{x}(t_k)\}\mathcal{H}_k\{\mathbf{x}(t_{k-1})\}^{\top}|\mathbf{y}_{1:K}, \hat{\mathbf{\psi}}^{(j-1)}\right],
\end{aligned}
$$

where the conditional expectations are computed in the E step via the EnKS.

# Appendix C

## Innovation likelihood

The innovation likelihood function is given by

$$ l(\mathbf{x}, \hat{\mathbf{\psi}}) = \prod_{k=1}^{K} \exp \left[ -\frac{1}{2} \mathbf{d}(t_k)^{\top} \left\{ \mathbf{P}_{yy}^{\mathrm{f}}(t_k) + \mathbf{R}(t_k) \right\}^{-1} \mathbf{d}(t_k) \right] \times (2\pi)^{-p/2} \left[ \det \left\{ \mathbf{P}_{yy}^{\mathrm{f}}(t_k) + \mathbf{R}(t_k) \right\} \right]^{-1/2}, $$
with the covariance matrix $P_{yy}^f(t_k)$ and the innovation vector $d(t_k)$ given in Eqs (12) and (13) respectively.

## References

Anderson JL. 2007. An adaptive covariance inflation error correction algorithm for ensemble filters. *Tellus* 59A: 210-224.

Annan JD, Hargreaves JC. 2004. Efficient parameter estimation for a highly chaotic system. *Tellus* 56A: 520-526.

Annan JD, Hargreaves JC. 2007. Efficient estimation and ensemble generation in climate modelling. *Phil. Trans. R. Soc. A: Math. Phys. Eng. Sci.* 365: 2077-2088.

Bougeault P, Jansa Clar A, Benech B, Carissimo B, Pelon J, Richard E. 1990. Momentum budget over the Pyrénées: The PYREX experiment. *Bull. Am. Meteorol. Soc.* 71: 806-818.

Burgers G, Van Leeuwen PJ, Evensen G. 1998. Analysis scheme in the ensemble Kalman filter. *Mon. Weather Rev.* 126: 1719-1724.

Cappé O, Moulines E, Rydøn T. 2005. *Inference in Hidden Markov Models*. Springer Science+Business Media: New York, NY.

Dee DP, Da Silva AM. 1998. Data assimilation in the presence of forecast bias. *Q. J. R. Meteorol. Soc.* 124: 269-295.

Dempster AP, Laird NM, Rubin DB. 1977. Maximum likelihood from incomplete data via the EM algorithm. *J. R. Stat. Soc. Ser. B* 39: 1-38.

Evensen G. 2009. *Data Assimilation: The Ensemble Kalman Filter*. Springer-Verlag: Berlin.

Evensen G, Van Leeuwen PJ. 2000. An ensemble Kalman smoother for nonlinear dynamics. *Mon. Weather Rev.* 128: 1852-1867.

Hertzog A, Alexander MJ, Plougonven R. 2012. On the intermittency of gravity wave momentum flux in the stratosphere. *J. Atmos. Sci.* 69: 3433-3448.

Hourdin F, Musat I, Bony S, Braconnot P, Codron F, Dufresne J-L, Fairhead L, Filiberti M-A, Friedlingstein P, Grandpeix J-Y, Krinner G, Levan P, Li Z-X, Lott F. 2006. The LMDZ4 general circulation model: Climate performance and sensitivity to parametrized physics with emphasis on tropical convection. *Clim. Dyn.* 27: 787-813.

Hu XM, Zhang F, Nielsen-Gammon JW. 2010. Ensemble-based simultaneous state and parameter estimation for treatment of mesoscale model error: A real-data study. *Geophys. Res. Lett.* 37: L08802, doi: 10.1029/2010GL043017.

Ide K, Courtier P, Ghil M, Lorenc AC. 1997. Unified notation for data assimilation: Operational, sequential and variational. *J. Meteorol. Soc. Jpn.* 75: 181-189.

Jackson C, Sen MK, Stoffa PL. 2004. An efficient stochastic Bayesian approach to optimal parameter and uncertainty estimation for climate model predictions. *J. Clim.* 17: 2828-2841.

Li F, Austin J, Wilson J. 2008. The strength of the Brewer-Dobson circulation in a changing climate: Coupled chemistry-climate model simulations. *J. Clim.* 21: 40-57.

Li H, Kalnay E, Miyoshi T. 2009. Simultaneous estimation of covariance inflation and observation errors within an ensemble Kalman filter. *Q. J. R. Meteorol. Soc.* 135: 523-533.

Liang X, Zheng X, Zhang S, Wu G, Dai Y, Li Y. 2011. Maximum likelihood estimation of inflation factors on error covariance matrices for ensemble Kalman filter assimilation. *Q. J. R. Meteorol. Soc.* 138: 263-273.

Lott F. 1995. Comparison between the orographic response of the ECMWF model and the PYREX 1990 data. *Q. J. R. Meteorol. Soc.* 121: 1323-1348.

Lott F. 1999. Alleviation of stationary biases in a GCM through a mountain drag parameterization scheme and a simple representation of mountain lift forces. *Mon. Weather Rev.* 127: 788-801.

Lott F, Miller MJ. 1997. A new subgrid-scale orographic drag parametrization: Its formulation and testing. *Q. J. R. Meteorol. Soc.* 123: 101-127.

Lott F, Fairhead L, Hourdin F, Levan P. 2005. The stratospheric version of LMDz: Dynamical climatologies, Arctic oscillation, and impact on the surface climate. *Clim. Dyn.* 25: 851-868.

McLandress C, Shepherd TG. 2009. Simulated anthropogenic changes in the Brewer-Dobson circulation, including its extension to high latitudes. *J. Clim.* 22: 1513-1540.

Miyoshi T. 2011. The Gaussian approach to adaptive covariance inflation and its implementation with the local ensemble transform Kalman filter. *Mon. Weather Rev.* 139: 1519-1535.

Miyoshi T, Kalnay E, Li H. 2012. Estimating and including observation-error correlations in data assimilation. *Inverse Prob. Sci. Eng.* 21: 387-398, doi: 10.1080/17415977.2012.712527.

NOAA. 2001. NOAA ETOPO2 Dataset, 2-minute Gridded Global Relief Data. National Geophysical Data Center: Boulder, CO. http://www.ngdc.noaa.gov/mgg/fliers/01mgg04.html (accessed 13 March 2014).

Palmer TN, Shutts GJ, Swinbank R. 1986. Alleviation of a systematic westerly bias in general circulation and numerical weather prediction models through an orographic gravity wave drag parametrization. *Q. J. R. Meteorol. Soc.* 112: 1001-1039.

Pham DT. 2001. Stochastic methods for sequential data assimilation in strongly nonlinear systems. *Mon. Weather Rev.* 129: 1194-1207.

Posselt DJ, Bishop CH. 2012. Nonlinear parameter estimation: Comparison of an ensemble Kalman smoother with a Markov chain Monte Carlo algorithm. *Mon. Weather Rev.* 140: 1957-1974.

Pulido M, Thuburn J. 2005. Gravity wave drag estimation from global analyses using variational data assimilation principles. I: Theory and implementation. *Q. J. R. Meteorol. Soc.* 131: 1821-1840.

Pulido M, Thuburn J. 2008. The seasonal cycle of gravity wave drag in the middle atmosphere. *J. Clim.* 21: 4664-4679.

Pulido M, Polavarapu S, Shepherd TG, Thuburn J. 2012. Estimation of optimal gravity wave parameters for climate models using data assimilation. *Q. J. R. Meteorol. Soc.* 138: 298-309.

Ruiz JJ, Pulido M, Miyoshi T. 2013. Estimating model parameters with ensemble-based data assimilation: A review. *J. Meteorol. Soc. Jpn.* 91: 79-99.

Schirber S, Klocke D, Pincus R, Quaas J, Anderson JL. 2013. Parameter estimation using data assimilation in an atmospheric general circulation model: From a perfect toward the real world. *J. Adv. Model. Earth Syst.* 5: 58-70.

Severijns CA, Hazeleger W. 2005. Optimizing parameters in an atmospheric general circulation model. *J. Clim.* 18: 3527-3535.

Sigmond M, Scinocca JF, Kushner PJ. 2008. Impact of the stratosphere on tropospheric climate change. *Geophys. Res. Lett.* 35: L12706, doi: 10.1029/2008GL033573.

Stainforth DA, Aina T, Christensen C, Collins M, Faull N, Frame DJ, Kettleborough JA, Knight S, Martin A, Murphy JM, Piani C, Sexton D, Smith LA, Spicer RA, Thorpe AJ, Allen MR. 2005. Uncertainty in predictions of the climate response to rising levels of greenhouse gases. *Nature* 433: 403-406.

Tandeo P, Ailliot P, Autret E. 2011. Linear Gaussian state-space model with irregular sampling: Application to sea surface temperature. *Stoch. Environ. Res. Risk Assess.* 25: 793-804.

Wang X, Bishop CH. 2003. A comparison of breeding and ensemble transform Kalman filter ensemble forecast schemes. *J. Atmos. Sci.* 60: 1140-1158.

Wu G, Zheng X, Wang L, Zhang S, Liang X, Li Y. 2013. A new structure for error covariance matrices and their adaptive estimation in EnKF assimilation. *Q. J. R. Meteorol. Soc.* 139: 795-804, doi: 10.1002/qj.2000.

Yang X, Delsole T. 2009. Using the ensemble Kalman filter to estimate multiplicative model parameters. *Tellus* 61A: 601-609.
## 4.5 Tandeo, Chapron, Ba, Autret et Fablet (2014) [TGRS]

**Contexte** Ce travail correspond à mon travail de recherche lorsque j’étais postdoctorant à Télécom Bretagne. Cette étude marqua le retour de mes recherches vers l’océanographie spatiale, en collaboration avec l’IFREMER et le LOPS. Plus précisément, c’est à ce moment là que j’ai découvert la synergie satellitaire, un champ de la télédétection s’intéressant à la corrélation entre variables issues de différents capteurs. Ce travail a été publié dans IEEE TGRS, un des journaux leader en télédétection spatiale. Cette étude fut à la base de la construction de l’école d’été ORSS (voir la section 1.3.3 pour plus de détails).

**Résumé** Les mesures satellitaires de la hauteur de la surface de la mer (SSH) et de la température de surface de la mer (SST), fournissent une grande quantité d’informations sur la circulation océanique, en particulier sur la dynamique des océans à méso-échelle, qui peut impliquer de fortes relations spatio-temporelles entre les champs de SSH et de SST. Dans un cadre purement basé sur les données, nous étudions dans quelle mesure la dynamique des océans à méso-échelle peut être décomposée en un mélange de modes dynamiques, caractérisés par différentes régressions locales entre les champs de SSH et de SST. Formellement, nous développons un nouveau modèle de régression basé sur des variables latentes, afin d’identifier des modes dynamiques cachés, à partir de séries d’observations conjointes de la SSH et de la SST. Appliqué à la région dynamique des Aiguilles, nous démontrons et discutons la pertinence géophysique du modèle de mélange, ainsi que sa capacité à réaliser une segmentation spatio-temporelle de la dynamique de la couche supérieure de l’océan.
# Segmentation of Mesoscale Ocean Surface Dynamics
Using Satellite SST and SSH Observations

Pierre Tandeo, Bertrand Chapron, Sileye Ba, Emmanuelle Autret, and Ronan Fablet

***Abstract***—Multisatellite measurements of altimeter-derived sea surface height (SSH) and sea surface temperature (SST) provide a wealth of information about ocean circulation, particularly mesoscale ocean dynamics which may involve strong spatiotemporal relationships between SSH and SST fields. Within an observation-driven framework, we investigate the extent to which mesoscale ocean dynamics may be decomposed into a mixture of dynamical modes, characterized by different local regressions between SSH and SST fields. Formally, we develop a novel latent class regression model to identify dynamical modes from joint SSH and SST observation series. Applied to the highly dynamical Agulhas region, we demonstrate and discuss the geophysical relevance of the proposed mixture model to achieve a spatiotemporal segmentation of the upper ocean dynamics.

***Index Terms***—Altimetry, remote sensing, sea surface, statistics, temperature.

## I. INTRODUCTION

IN THE last two decades, multisatellite measurements of altimeter-derived sea surface height (SSH) and multisensor measurements of sea surface temperature (SST) have provided a wealth of information about ocean circulation and atmosphere-ocean interactions. As a depth-integrated quantity dependent upon the density structure of the water column, altimeter SSH estimations capture mesoscale structures, horizontal scales of 50 km to few hundred kilometers, and allow for the retrieval of surface currents using the geostrophy balance. This emerging and rich mesoscale circulation further stirs the large-scale SST fields. Accordingly, our picture of upper ocean dynamics has considerably evolved toward a complex system characterized by strong interactions, whose spatiotemporal variability extends over a wide range of scales. Furthermore, several studies (cf. [13], [17], [19], [20], and [21]) rationalize and demonstrate that fields of SST can become an active tracer coupled to the dynamics, leading to strong correlations with SSH fields.

Such a framework can possibly guide the investigation and implementation of improved statistical means to optimally combine existing multialtimeter SSH measurements with other satellite medium- to high-resolution observations (e.g., microwave SST and salinity and scatterometer winds), augmented by the growing available *in situ* data (e.g., [1], [8], and [26]). Theoretically, the upper ocean turbulence for the horizontal scales between 50 km to few hundred kilometers is still consistent with the geostrophy turbulence theory. Under this assumption, the upper ocean dynamics may be simply predicted from surface density horizontal variations possibly dominated by SST variations. For such a case, a linear transfer function shall be identified between SSH and SST fields to also possibly lead to the estimation of the subsurface flow (e.g., [16] and [20]). This linear transfer function does not involve temporal differencing as in the maximum cross-correlation technique or alternate strategies (e.g., [4], [7], and [23]). Note that other recent papers use nonlinear transfer functions to relate SST and SSH fields (e.g., [12]).

This strongly advocates for observation-driven studies to explore and characterize the local relationships between SST, SSH, and the derived surface currents from satellite-based routine observations. However, as illustrated in Fig. 1, a simple linear transfer function cannot be expected to solely govern the whole mesoscale dynamics in a particular ocean region. As revealed, an overall spatial correlation exists, but for instance, relationships between SST gradients and altimetry-derived surface currents may spatially differ. In the warmer SST frontal zone, SST gradients correspond to large surface currents (top of the image). In the colder frontal area, large SST gradients do not reflect in large surface currents (bottom of the image). Moreover, the clearly detected eddy (top left of the image) is associated with weak SST gradients.

Within an observation-driven framework, one may consider joint principal component analysis/empirical orthogonal function (EOF) procedures to decompose the relationships between SST and SSH fields, as widely used in ocean sensing applications (cf. [6] and [22]). Such EOF-based schemes would, however, resort to a single linear and global model. As such, this model could not address the spatial nonstationarity of the SST-SSH relationships illustrated in Fig. 1. By contrast, we here consider local linear transfer functions. We assume that, locally, upper ocean dynamics may be analyzed according to a finite mixture model, where each component of the mixture is characterized by a local SST-SSH linear transfer function. This mixture-based representation relates to a nonlinear model with assumptions in accordance with the observations made earlier between SST and surface currents. In this paper, we propose to investigate such a model to accomplish the following: 1) to

Manuscript received January 3, 2013; revised April 5, 2013 and July 18, 2013; accepted August 23, 2013. This work was supported in part by Grant F120, the LabexMER, Agence Nationale de la Recherche RedHots, and Institut Mines-Telecom.

P. Tandeo and R. Fablet are with Telecom Bretagne, 29280 Plouzané, France.

B. Chapron and E. Autret are with IFREMER, 29280 Plouzané, France.

S. Ba is with the RN3D Innovation Lab, 13003 Marseille, France.

This paper has supplementary downloadable material available at http://ieeexplore.ieee.org, provided by the authors. It shows animations of the time series of daily maps in the Agulhas current over 2004 shown in Fig. 4.

Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org.

Digital Object Identifier 10.1109/TGRS.2013.2280494
<figure><img src="image_72.png" /><figcaption>Fig. 1. Surface currents derived from altimeter SSHs and (a) microwave SSTs and (b) the associated temperature gradient norms within the Agulhas return current on January 1, 2004.</figcaption></figure>

develop a probabilistic learning-based setting for the inference of such mixture models and the spatiotemporal segmentation of the identified dynamical modes (i.e., the different components of the mixture model) and 2) to evaluate the extent to which such mixture models are geophysically relevant to characterize the upper ocean dynamics over active ocean regions.

Hereafter, we consider the Agulhas region, and this paper is organized as follows. Section II presents the remote sensing data and describes our probabilistic learning-based model. In Section III, the application to satellite observations is evaluated. We further discuss and summarize the key results of our investigations in Section IV.

## II. DATA AND METHODOLOGY

### A. Patch-Based Approach

As mentioned earlier, recent theoretical and numerical experiments have stressed that upper ocean dynamics may be characterized by couplings between SSH and SST according to the following relationship in the Fourier domain (cf. [14]):

$$ \mathcal{F}_H(\tilde{\text{SSH}}) = -\gamma |\mathbf{k}|^{-\alpha} \mathcal{F}_T(\tilde{\text{SST}}) \quad (1) $$

<figure><img src="image_73.png" /><figcaption>Fig. 2. Sketch of SST patches (in degrees represented in false colors), denoted by $\mathbf{X}$, and the corresponding SSHs (in meters represented by dots) and surface currents (in meters per second represented by quivers) denoted by $\mathbf{Y}$ at the central location $s_i$ and $s_j$ at time $t_i$.</figcaption></figure>

where $k$ is the horizontal wavenumber vector and $\mathcal{F}_H$ and $\mathcal{F}_T$ are linear filters of SSH and SST, respectively. The $\alpha$ parameter sets up the effective coupling between surface fields. For $\alpha = 1$, (1) resorts to the surface quasi-geostrophy (SQG) model. In [16], $\mathcal{F}_H$ and $\mathcal{F}_T$ were bandpass filters between 80 and 300 km. As $\alpha$ increases, the smoothing increases, and couplings decrease. For $\alpha = 2$, the filtered SST would trace the vorticity. Formally, (1) states that surface currents can be regarded as spatial derivatives of a filtered version of the SST field. The parameter $\gamma$ relates to a normalization constraint. In general, parameters $\gamma$ and $\alpha$, as well as the definition of the filters $\mathcal{F}_H$ and $\mathcal{F}_T$, may spatially vary such that a single linear transfer function as in (1) is unlikely to apply globally as illustrated in Fig. 1.

These considerations led us to hypothesize that zonal and meridional surface currents noted (U, V) and SSH can still locally relate to SST derivatives, but according to a finite set of $K$ linear transfer functions, hereafter referred to as $K$ dynamical modes. Formally, this is stated in the Fourier domain as

$$ (\tilde{\text{SSH}}, \tilde{\text{U}}, \tilde{\text{V}}) = \mathcal{H}_k(\tilde{\text{SST}}) \quad (2) $$

where $\mathcal{H}_k$ characterizes the $k$th dynamical mode, which locally relates SSH and SST fields through linear filter $\mathcal{H}_k$. In this paper, we do not consider any bandpass filters $\mathcal{F}_H$ and $\mathcal{F}_T$. Using a matrix formulation, (2) is rewritten in the real domain as a patch-based linear regression

$$ \mathbf{Y}(s_i, t_i) = \mathbf{H}_k(\mathbf{X}(s_i, t_i)) \quad (3) $$

where $\mathbf{Y}(s_i, t_i)$ encodes the local SSH variability through a 3-D vector formed by the SSH value and the surface current (U, V) at spatiotemporal location $(s_i, t_i)$ and $\mathbf{X}(s_i, t_i)$ is the vectorized version of the local SST patch centered in $s_i$ at time $t_i$ (cf. Fig. 2). It may be noted that we encoded local SSH variations at spatiotemporal location $(s_i, t_i)$ through the surface currents which are computed as the spatial derivatives of the SSH field. As such, it constrains the method to account for spatial regularity. The linear operator associated with dynamical mode $k$ is
corresponding to the local version of $\mathcal{H}_k$ in (2). It corresponds to three vectorized versions of spatial convolution matrices. Here, $p$ defines the size of the local SST neighborhood around $s_i$ and is set according to the Rossby radius of the study region, i.e., the mean size of the mesoscale ocean structures like eddies. For the Great Agulhas current region, we set it up to 200 km, i.e., $p = 81$ for the spatial resolution of the considered data.

## B. Remote Sensing Data

As SSH and surface geostrophic current (U, V) data, we use the daily delayed time Maps of Absolute Dynamic Topography (MADT) produced by Collecte Localisation Satellites, available online at http://www.aviso.oceanobs.com/. This information combines the signal of several altimeters onto a 1/3° Mercator projection grid. We use the 2004 data since four altimeters were available (Jason-1, Envisat or European Remote-Sensing-2, Topex/Poseidon, and Geosat Follow-On). As SST data, we use optimally interpolated microwave SSTs provided by Remote Sensing System (RSS), available online at http://www.ssmi.com/. It combines the signal of three microwave radiometers (Tropical rainfall measuring missions Microwave Imager, Advanced Microwave Scanning Radiometer Earth observing system, and WindSAT) which are robust to the presence of clouds. The spatial resolution is 1/4° × 1/4°, and the temporal resolution is the same as the MADT data, i.e., daily. We bilinearly interpolate the MADT data onto the SST grid. We focus on the Agulhas region between longitudes 5° E to 65° E and latitudes 30° S to 48° S.

Given the joint series of satellite observations, we extract SST patches (denoted by $\mathbf{X}$) and the associated SSH with surface current (U, V) at the center of the patches (denoted by $\mathbf{Y}$). Overall, the processed data set is composed of ~5 × 10⁶ pairs of vectors $\mathbf{X}$ and $\mathbf{Y}$. To infer the parameters of the considered mixture model, i.e., the parameters of each dynamical mode in (3), we first build a training data set as a random sample of $n = 10^5$ elements (for a given day, we use about 2% of the data to fit the model). In a second step, we apply the inferred mixture model to the entire processed data set to extract the spatiotemporal of the different dynamical modes.

## C. Latent Class Regression Model

Our objective is to identify $K$ different dynamical modes (latent variable $Z$) from a joint set of SST patches ($p$-dimensional vector $\mathbf{X}$) and SSH with zonal and meridional surface currents (3-D vector $\mathbf{Y}$). In this paper, we assume that the conditional probability of $\mathbf{Y}$ given $\mathbf{X}$ and the dynamical mode $Z = k$ is given by

$$p(\mathbf{Y}|\mathbf{X}, Z = k) \propto \mathcal{N}_k(\mathbf{Y}; \mathbf{X}\beta_k, \Sigma_k) \quad (4)$$

where $\mathcal{N}_k$ represents a multivariate Gaussian probability density function evaluated in $\mathbf{Y}$ with mean $\mathbf{X}\beta_k$ and covariance $\Sigma_k$. Hence, the conditional probability of $\mathbf{Y}|\mathbf{X}$ resorts to a mixture of normal distributions

$$p(\mathbf{Y}|\mathbf{X}, \boldsymbol{\theta}) = \sum_{k=1}^{K} \lambda_k \mathcal{N}_k(\mathbf{Y}; \mathbf{X}\beta_k, \Sigma_k) \quad (5)$$

where $\lambda_k$ is the prior probability of mode $k$. To simplify the notations, we store the overall parameters of (5) in $\boldsymbol{\theta} = (\lambda_1, \dots, \lambda_K, \beta_1, \dots, \beta_K, \Sigma_1, \dots, \Sigma_K)$. In the literature, this model is referred to as a “latent class regression” or “cluster-wise regression” (cf. [10]). By construction, it imposes that $0 \le \lambda_k \le 1$, $\sum_{k=1}^K \lambda_k = 1$, and $\Sigma_k$ is positive defined. The maximum likelihood estimation procedure for model parameters $\boldsymbol{\theta}$ is given hereinafter.

## D. Model Learning

To learn model parameters $\boldsymbol{\theta}$ in (5), we resort to a maximum likelihood criterion and use an iterative expectation-maximization (EM) procedure (cf. [9]). It relies on the maximization of the log-likelihood given by

$$\mathcal{L}(\boldsymbol{\theta}) = \sum_{i=1}^{n} \log(p(\mathbf{Y}(s_i, t_i)|\mathbf{X}(s_i, t_i), \boldsymbol{\theta})) \quad (6)$$

where $n$ is the number of observations of the training data set described in Section II-B. From a given initialization, the EM procedure iterates an expectation step (E-step) and a maximization step (M-step). At a given iteration, using the Bayes theorem, the E-step resorts to the computation of the posterior probabilities of the latent variable $Z$ for each spatiotemporal location $(s_i, t_i)$ given the current parameter estimate $\hat{\boldsymbol{\theta}}$

$$\hat{\pi}_k(s_i, t_i) = \frac{\hat{\lambda}_k \mathcal{N}_k(\mathbf{Y}(s_i, t_i); \mathbf{X}(s_i, t_i)\hat{\beta}_k, \hat{\Sigma}_k)}{p(\mathbf{Y}(s_i, t_i)|\mathbf{X}(s_i, t_i), \hat{\boldsymbol{\theta}})}, \quad \forall k. \quad (7)$$

The M-step then minimizes the expectation of the log-likelihood conditionally to the current parameter estimate $\hat{\boldsymbol{\theta}}$. This leads to the update of the prior probabilities of the latent variable $Z$ as

$$\hat{\lambda}_k = \frac{\sum_{i=1}^{n} \hat{\pi}_k(s_i, t_i)}{n}, \quad \forall k. \quad (8)$$

The updated regression parameters $\hat{\beta}_k$, $\forall k$ are derived by fitting $K$ distinct linear regressions using a weighted least square criterion on the $n$ observations where the weights are given by the posterior probabilities given in (7) as in [25]. Then, we maximize $\mathcal{L}$ with respect to $\Sigma_k$ and obtain

$$\hat{\Sigma}_k = \frac{\sum_{i=1}^{n} \hat{\pi}_k(s_i, t_i) \epsilon_k(s_i, t_i)^{\top} \epsilon_k(s_i, t_i)}{\sum_{i=1}^{n} \hat{\pi}_k(s_i, t_i)}, \quad \forall k \quad (9)$$

where $\epsilon_k(s_i, t_i) = \mathbf{Y}(s_i, t_i) - \mathbf{X}(s_i, t_i)\hat{\beta}_k$ and $^{\top}$ represents the transpose. The algorithm iterates the E-step and M-step until a negligible increase of the log-likelihood $\mathcal{L}$ which is strictly growing.

A key aspect of the latent class regression is the choice of $K$, the number of clusters. Different statistical criteria state the selection of parameter $K$ as a tradeoff between the likelihood and the complexity of the model (cf. [18]) such as the Akaike information criterion (AIC). However, the optimization of these criteria makes no effort to distinguish the error explained by the regression fit and the error explained by the clustering process. In practice, there is an actual potential for overfitting
with the latent class regression model (cf. [5]). Then, for a given number of clusters K, we suggest different evaluations of the EM algorithm to reach a greater consistency in the estimation of model parameters (cf. [15]). The idea is to use random values $\hat{\pi}_k$ as initialization values of the EM procedure and select parameter estimates corresponding to the greatest likelihood (see [3] for more details). In our case, the inference based on AIC would lead to K values between 4 and 9. We performed a complementary qualitative analysis, and the setting K = 4 resulted in a good tradeoff between the geophysical interpretation of the model and regression error statistics (i.e., the maximization of the likelihood).

### E. Spatiotemporal Segmentation of Dynamical Modes and SSH/Current Predictions

We exploit the inferred mixture model with parameter $\hat{\theta}$ to perform a spatiotemporal segmentation of the underlying dynamical modes. More precisely, for any spatial location s and time t within the Agulhas region over the year 2004 (cf. processed data in Section II-B), we use (7) to evaluate the posterior probability $\hat{\pi}_k(s, t)$ for the K = 4 dynamical modes. Then, the pixels will be assigned to the most likely dynamical mode. One can also estimate for each time t the relative spatial occurrence of each dynamical mode using (8). Moreover, using (3), the estimation of the SSH and surface current at the spatial location s and time t from the associated SST patch follows from the fuzzy regression

$$ \hat{\mathbf{Y}}(s, t) = \sum_{k=1}^{K} \hat{\pi}_k(s, t) \mathbf{X}(s, t) \hat{\beta}_k \quad (10) $$

where the $\hat{\beta}_k$ parameterize the linear operator $\mathbf{H}_k$. They correspond to the $p \times 3$ regression coefficient matrices between the patch of SSTs ($\mathbf{X}$ in degrees) and the central value of SSH along with zonal U and meridional V components of the surface currents ($\mathbf{Y}$ in meters and meters per second). All the locations of the Agulhas current share the same K matrices $\hat{\beta}_k$ which are not changing temporally nor spatially. It may be outlined that no additional constraint is set on the posterior probabilities $\hat{\pi}_k(s, t)$, which could reveal space-time variations of the distribution of the dynamical modes (including seasonal variations).

### III. CHARACTERIZATION OF OCEAN SURFACE DYNAMICS

We first report the temporal evolution of the relative spatial occurrence of the K = 4 dynamical modes (cf. Fig. 3). The dynamical modes involve clear seasonal cycles. Dynamical modes 1 (red circles) and 4 (blue diamonds) depict similar temporal variations, completely out of phase with dynamical modes 2 (green crosses) and 3 (cyan squares). To study the spatial distribution of these dynamical modes, we focus on two dates corresponding to the maximal and minimal values of the seasonal cycle, namely, March 1 and September 1, 2004. They correspond to the maximum or minimum values of the $\hat{\lambda}_k$.

For these two dates, from the maps of posterior probabilities $\hat{\pi}_k(s, t)$, we determine the segmentation maps of each

<figure><img src="image_74.png" /><figcaption>Fig. 3. Time series of the relative proportion of locations associated with the four dynamical modes within the Agulhas region for the year 2004. In the subsequent, red circles, green crosses, cyan squares, and blue diamonds respectively distinguish the first, second, third, and fourth dynamical modes. The straight lines correspond to a six-month time separation: March 1 and September 1, 2004. The corresponding spatial maps of the posterior probabilities for each dynamical mode are given in Fig. 4.</figcaption></figure>

dynamical mode as illustrated in Fig. 4. The animations of the time series of these daily maps in the Agulhas current over 2004 are available as supplementary material or online at http://tandeo.wordpress.com/communications/articles/. A qualitative analysis of these maps highlights a clear spatiotemporal segmentation of the different dynamical modes that can be interpreted from a geophysical perspective in terms of different geophysical processes. We also report for each dynamical mode the observed distributions of current, height, and temperature values (cf. Fig. 5). The first dynamical mode (red circles) characterizes very strong current magnitude and warm waters. It is primarily associated with the main Agulhas current that flows down the east coast of Africa through the Agulhas ridge. This mode also involves mesoscale eddies, the so-called warm core Agulhas rings with strong surface currents, low temperature gradients, and middle-range SSH values around 0.5 m. The latter seems to be a discriminative feature of this first mode. The second dynamical mode (green crosses) mainly relates to the eastward Agulhas return current that hits a part of the South Atlantic current. It creates a subtropical front varying from 36° S to 44° S with strong eastward currents, middle-range SST gradients, and large SSH values (about 1 m). This front was clearly observed in the upper part of Fig. 1(b). The third (cyan squares) and fourth (blue diamonds) dynamical modes correspond to weaker surface currents. The third one is characterized by midtemperatures and westward currents, whereas the fourth one involves colder temperatures and eastward currents. Let us stress that the third dynamical mode involves large SST gradients but weak surface currents as identified in the lower part of Fig. 1(b). In this mode, the SST can be clearly identified as a passive tracer of the surface upper ocean dynamics.

To characterize more precisely the inferred model, we plot the regression line and the 95% confidence region from the estimated parameters $\hat{\beta}_k$ and $\hat{\Sigma}_k$ of each dynamical mode (cf. Fig. 6). These results clearly stress the relevance of a mixture model, compared to a single linear model (represented as the fine black line). Similar slopes are observed for the second and third dynamical modes as well as the first and fourth ones. Stronger differences among the dynamical modes are outlined regarding the associated regression error variances. We can notice that the first mode involves the greater variance and the fourth one involves the lower one.
<figure><img src="image_75.png" /><figcaption>Fig. 4. Maps of the posterior probabilities given in (7) of the dynamical modes given the SST and SSH fields on (a) March 1 and (b) September 1, 2004 within the Agulhas current. We use a four-class latent regression model fitted from the whole year of 2004 (see text for details). For a given location and time, the sum of the four probabilities is equal to 1. The animations of the daily maps are available as supplementary material.</figcaption></figure>

We further investigate the geophysical consistency of the identified dynamical modes from the comparison between the true MADT SSH and surface currents in Fig. 7(a) and the prediction of the latent class regression model in Fig. 7(b) given by (10). Overall, as reported in Table I, a good agreement is obtained with a global correlation coefficient of 0.72 for the surface currents and 0.96 for the SSH; it can locally be very large as illustrated in the right column of Fig. 7 (corresponding to the zone depicted in Fig. 1). This zone involves the four dynamical modes. The mixture model enables us to retrieve both the large warm eddy associated with weak SST gradients (top left of the zone and first dynamical mode) and the relatively large surface currents along the large warmer SST gradients (upper part and second mode) as well as the rather weak currents along the large but colder SST front (lower part and third mode). For comparison purposes, we also plot the results issued from a single linear transfer function in Fig. 7(c). This

model clearly underestimates the surface currents within the warm eddy (top left) and overestimates the currents of the colder frontal zone (lower part), which stresses the requirement for considering a mixture model.

To further characterize each dynamical mode, we report in Table II the correlation and root mean square error (RMSE) statistics computed within the associated spatiotemporal domain, i.e., the domain comprising all spatiotemporal locations assigned to dynamical mode *k* according to posterior probabilities $\hat{\pi}_k(s, t)$ computed in (7). We compare the latent class regression model with respect to both the true MADT data and an SQG-like hypothesis, i.e., (1) with $\alpha = 1$, within the space-time region associated with each dynamical mode. This analysis clearly discriminates the second and fourth modes from the first and third modes. The linear transfer functions of the second and fourth dynamical modes capture a large part of the variability of the true SSH data. These first two modes
<figure><img src="image_76.png" /></figure>

<figure><img src="image_77.png" /></figure>

<figure><img src="image_78.png" /></figure>

<figure><img src="image_79.png" /><figcaption>Fig. 5. Distribution of the (a) SST, (b) SSH, and surface current (U, V) (c) norm and (d) direction. The results are given for each dynamical mode within the Agulhas current for the whole year of 2004.</figcaption></figure>

have also a good consistency with the SQG hypothesis, with correlation coefficients of 0.63 and 0.68, respectively. These results suggest that the SST might be regarded as an active tracer of the surface dynamics in the associated regions. By contrast, the SQG hypothesis poorly fits to the first and third dynamical modes, with correlation coefficients of 0.33 and 0.27, respectively. These two dynamical modes also involve a slightly lower predictability of the linear transfer functions to retrieve the SSH and surface currents.

To summarize, the four dynamical modes correspond to different physical parameter values in (1). The second and fourth modes seem to correspond to $\alpha \approx 1$ (close to the SQG model), whereas the first and fourth modes appear to be $\alpha > 1$ (coupling of SST and SSH at large scales). The factor $\gamma$ relates to the amplitude of the surface geostrophic currents: Large values relate to strong currents as retrieved for modes 1 and 2, and low values relate to strong currents as observed for modes 3 and 4. Overall, these results are consistent with the previous work reported in [16] and [27]. In particular,

<figure><img src="image_80.png" /></figure>

<figure><img src="image_81.png" /><figcaption>Fig. 6. SSH as a function of SST. For each dynamical mode, we give the regression line and the 95% confidence envelope. The four regressions are highly significant (p-values < $10^{-3}$). The fine black line is a benchmark curve corresponding to the global linear regression (with a one-mode model).</figcaption></figure>

Isern-Fontanet et al. explored the SQG hypothesis from a phase-correlation analysis between SSH and SST fields, while Xu et al. explored the coupling assumption from a spectral analysis of SSH fields. These authors concluded that SQG-like dynamics would mainly occur near the edge of the large current system. Compared to these analysis, our contribution is twofold: the quantitative characterization of the extent to which the SQG dynamics applies through correlation statistics as well as the actual space-time tracking of the regions associated with SQG-like and non-SQG dynamical modes.

## IV. CONCLUSION AND PERSPECTIVES

In this paper, we propose an observation-driven framework to identify, characterize, and track ocean surface dynamical modes. We rely on a latent class regression model, where the dynamical modes are characterized by a local linear transfer function between SST, SSH, and surface current (U, V), in agreement with the theoretical assumption given in (1). This probabilistic approach locally relates the distribution of the SSH and sea surface currents conditionally to the SST via a nonlinear model: a Gaussian mixture of linear transfer functions. The statistical parameters of the model are estimated using a maximum likelihood approach.

We applied the proposed methodology to the 2004 daily $1/4^\circ \times 1/4^\circ$ satellite SST and SSH image series. The reported results retrieved a relevant spatiotemporal decomposition of ocean surface dynamics in the Agulhas region according to four dynamical modes: 1) the main Agulhas current and warm core rings characterized by strong currents and hot temperatures where the SST is weakly correlated with SSH; 2) the return Agulhas current with lower temperatures and currents where the SST is an active tracer; 3) local front regions where strong SST gradients do not seem to affect the current velocities; and 4) a weaker dynamical mode where SST is strongly correlated to SSH.
<figure><img src="image_82.png" /><figcaption>Fig. 7. (a) True MADT data and the results for (b) the proposed latent class regression model using (10) and for (c) the linear model on January 1, 2004 within the Agulhas current. The left column corresponds to the SSH, and the right column corresponds to the SSH and surface currents (U, V) for the zone depicted in Fig. 1 and in the black box.</figcaption></figure>

Our study complements previous theoretical studies which showed that mesoscale upper ocean dynamics may be characterized by a linear coupling between SST and SSH (cf. [13], [17], [19], [20], and [21]). Following a fully observation-driven framework, the proposed latent regression model enabled us to identify different dynamical modes, including some SQG-like ones, and to track the space-time extension of each dynamical mode. The reported results clearly pointed out the requirement for considering a mixture model to decompose the space-time variabilities of the ocean surface dynamics. Regarding methodological aspects, it may be pointed out that EOF-based schemes (cf. [6] and [22]) could not reveal such nonstationary space-time variabilities. The joint EOF scheme typically decomposes a global linear mapping between the analyzed fields according to principal modes. In our case, such an EOF decomposition could be considered for each dynamical mode to further characterize the associated linear transfer function with respect to joint SST-SSH principal modes.

Regarding our future work, we will further investigate latent class regression models with additional regressors. Among others, it seems appealing to explore how time-lagged SST features and other geophysical fields such as wind speed, mixed-layer

<table><caption>TABLE I
CORRELATION AND RMSE STATISTICS WITHIN THE AGULHAS CURRENT FOR THE WHOLE YEAR OF 2004. THE LABELS "MADT," "LINEAR," AND "LATENT" REFER RESPECTIVELY TO THE TRUE MADT DATA, THE LINEAR MODEL, AND THE PROPOSED LATENT CLASS REGRESSION MODEL. ALL THE CORRELATIONS ARE STATISTICALLY SIGNIFICANT (P-VALUES < 10⁻³)</caption><thead><tr><th colspan="2">Correlation (RMSE)</th><th>SSH</th><th>(U,V)</th></tr></thead><tbody><tr><td>MADT</td><td>LINEAR</td><td>0.89 (0.22)</td><td>0.61 (0.19)</td></tr><tr><td>MADT</td><td>LATENT</td><td>0.96 (0.16)</td><td>0.72 (0.16)</td></tr></tbody></table>

depth, salinity, and chlorophyll-a concentration (cf. [24]) could improve the estimation of SSH and surface currents. We also plan to apply the proposed model to other strongly active ocean regions such as the Gulf Stream system. Our objective will be to determine shared and/or system-specific dynamical modes. Future work will also investigate more detailed physical interpretation of the identified dynamical modes, particularly in terms of spectral characteristics. For instance, it would be interesting to relate more precisely the physical parameters $\gamma$ and $\alpha$ of (1) to the different hidden dynamical modes extracted by our statistical approach. Whereas factor $\gamma$ seems to be well estimated in this paper, the spatial resolution of the satellite
<table><caption>TABLE II
CORRELATION AND RMSE STATISTICS FOR SPATIOTEMPORAL LOCATIONS ASSIGNED TO EACH DYNAMICAL MODE ACCORDING TO POSTERIOR PROBABILITIES COMPUTED IN (7) WITHIN THE AGULHAS CURRENT FOR THE WHOLE YEAR OF 2004. THE LABELS "MADT," "SQG," AND "LATIVE" REFER RESPECTIVELY TO THE TRUE MADT DATA, A SURFACE QUASI-GEOSTROPHIC HYPOTHESIS, AND THE PROPOSED LATENT CLASS REGRESSION MODEL. ALL THE CORRELATIONS ARE STATISTICALLY SIGNIFICANT (P-VALUES < 10⁻³)</caption><thead><tr><th>Correlation (RMSE)</th><th>Mode 1</th><th>Mode 2</th><th>Mode 3</th><th>Mode 4</th></tr></thead><tbody><tr><td>SSH<sub>MADT</sub></td><td>0.72 (0.26)</td><td>0.93 (0.16)</td><td>0.97 (0.13)</td><td>0.96 (0.12)</td></tr><tr><td>(U,V)<sub>MADT</sub></td><td>0.62 (0.33)</td><td>0.88 (0.15)</td><td>0.63 (0.08)</td><td>0.88 (0.06)</td></tr><tr><td>(U,V)<sub>SQG</sub></td><td>0.33 (0.41)</td><td>0.63 (0.28)</td><td>0.27 (0.25)</td><td>0.68 (0.18)</td></tr></tbody></table>

data (up to 50 km) does not permit to detect dynamical modes with $\alpha < 1$, corresponding to more local couplings between geophysical fields. Moreover, such an improved model shall then possibly address both of the following: 1) the higher resolution prediction of mesoscale ocean surface currents from SST spatiotemporal fields (and improve the nongeostrophic component estimation as in [7] and [23]) and 2) the extraction of new local and global descriptors of ocean surface dynamics from satellite sea surface observations (cf. [2] and [11]).

ACKNOWLEDGMENT

The authors would like to thank the Archiving, Validation and Interpretation of Satellite Oceanographic and the Remote Sensing System projects for respectively providing the altimeter-derived SSH and surface current and the microwave SST data. The authors would also like to thank the reviewers and the editor for their numerous constructive comments.

REFERENCES

[1] M. M. Ali, D. Swain, and R. A. Weller, "Estimation of ocean subsurface thermal structure from surface parameters: A neural network approach," *Geophys. Res. Lett.*, vol. 31, no. 20, pp. L20308-1-L20308-4, Oct. 2004.

[2] S. O. Ba, E. Autret, B. Chapron, and R. Fablet, "Statistical descriptors of ocean regimes from the geometric regularity of SST observations," *IEEE Geosci. Remote Sens. Lett.*, vol. 9, no. 5, pp. 851-855, Sep. 2012.

[3] C. Biernacki, G. Celeux, and G. Govaert, "Choosing starting values for the EM algorithm for getting the highest likelihood in multivariate Gaussian mixture models," *Comput. Statist. Data Anal.*, vol. 41, no. 3/4, pp. 561-575, Jan. 2003.

[4] M. M. Bowen, W. J. Emery, J. L. Wilkin, P. C. Tildesley, I. J. Barton, and R. Knewtson, "Extracting multiyear surface currents from sequential thermal imagery using the maximum cross-correlation technique," *J. Atmos. Ocean. Technol.*, vol. 19, no. 10, pp. 1665-1676, Oct. 2002.

[5] M. J. Brusco, J. D. Cradit, D. Steinley, and G. L. Fox, "Cautionary remarks on the use of clusterwise regression," *Multivariate Behavioral Res.*, vol. 43, no. 1, pp. 29-49, 2008.

[6] K. S. Casey and D. Adamec, "Sea surface temperature and sea surface height variability in the North Pacific Ocean from 1993 to 1999," *J. Geophys. Res.*, vol. 107, no. C8, pp. 14-1-14-12, 2002.

[7] W. Chen, "Surface velocity estimation from satellite imagery using displaced frame central difference equation," *IEEE Trans. Geosci. Remote Sens.*, vol. 50, no. 7, pp. 2791-2801, Jul. 2012.

[8] P. C. Chu, C. Fan, and W. T. Liu, "Determination of vertical thermal structure from sea surface temperature," *J. Atmos. Ocean. Technol.*, vol. 17, pp. 971-979, 2000.

[9] A. P. Dempster, N. M. Laird, and D. B. Rubin, "Maximum likelihood from incomplete data via the EM algorithm," *J. R. Stat. Soc. Series B (Methodological)*, vol. 39, no. 1, pp. 1-38, 1977.

[10] W. S. DeSarbo and W. L. Cron, "A maximum likelihood methodology for clusterwise linear regression," *J. Classification*, vol. 5, no. 2, pp. 249-282, Sep. 1988.

[11] R. Fablet, A. Chaigneau, and S. Bertrand, "Multiscale geometric deformations along planar curves: Application to satellite tracking and ocean observation data," *IEEE Trans. Geosci. Remote Sens.*, to be published.

[12] L. Gaultier, J. Verron, J. M. Brankart, O. Titaud, and P. Brasseur, "On the inversion of submesoscale tracer fields to estimate the surface ocean circulation," *J. Marine Syst.*, vol. 126, pp. 33-42, Oct. 2013.

[13] U. Hausmann and A. Czaja, "The observed signature of mesoscale eddies in sea surface temperature and the associated heat transport," *Deep Sea Res. Part I, Oceanogr. Res. Papers*, vol. 70, pp. 60-72, Dec. 2012.

[14] I. M. Held, R. T. Pierrehumbert, S. T. Garner, and K. L. Swanson, "Surface quasi-geostrophic dynamics," *J. Fluid Mech.*, vol. 282, pp. 1-20, 1995.

[15] C. Hennig, "Identifiability of models for clusterwise linear regression," *J. Classification*, vol. 17, pp. 273-296, Jul. 2000.

[16] J. Isern-Fontanet, B. Chapron, G. Lapeyre, and P. Klein, "Potential use of microwave sea surface temperatures for the estimation of ocean currents," *Geophys. Res. Lett.*, vol. 33, no. 24, pp. L24608-1-L24608-5, Dec. 2006.

[17] J. Isern-Fontanet, G. Lapeyre, P. Klein, B. Chapron, and M. W. Hecht, "Three-dimensional reconstruction of oceanic mesoscale currents from surface information," *J. Geophys. Res.*, vol. 113, no. C9, pp. C09005-1-C09005-17, Sep. 2008.

[18] J. B. Kadane and N. A. Lazar, "Methods and criteria for model selection," *J. Amer. Statist. Assoc.*, vol. 99, no. 465, pp. 279-290, Mar. 2004.

[19] P. Klein, J. Isern-Fontanet, G. Lapeyre, G. Roullet, E. Danioux, B. Chapron, S. Le Gentil, and H. Sasaki, "Diagnosis of vertical velocities in the upper ocean from high resolution sea surface height," *Geophys. Res. Lett.*, vol. 36, pp. L12603-1-L12603-5, 2009.

[20] H. LaCasce and A. Mahadevan, "Estimating subsurface horizontal and vertical velocities from sea-surface temperature," *J. Marine Res.*, vol. 64, no. 5, pp. 695-721, Sep. 2006.

[21] G. Lapeyre and P. Klein, "Dynamics of the upper oceanic layers in terms of surface quasigeostrophy theory," *J. Phys. Oceanogr.*, vol. 36, pp. 165-176, Feb. 2006.

[22] E. W. Leuliette and J. M. Whar, "Coupled pattern analysis of sea surface temperature and TOPEX/Poseidon sea surface height," *J. Phys. Oceanogr.*, vol. 29, pp. 599-611, 1999.

[23] J. Marcello, F. Eugenio, F. Marqués, A. Hernandez-Guerra, and A. Gasull, "Motion estimation techniques to automatically track oceanographic thermal structures in multisensor image sequences," *IEEE Trans. Geosci. Remote Sens.*, vol. 46, no. 9, pp. 2743-2762, Sep. 2008.

[24] M. Saraceno, C. Provost, and A. R. Piola, "On the relationship between satellite-retrieved surface temperature fronts and chlorophyll a in the western South Atlantic," *J. Geophys. Res.: Oceans (1978-2012)*, vol. 110, no. C11, pp. C11016-1-C11016-18, Nov. 2005.

[25] R. Schlittgen, "A weighted least-squares approach to clusterwise regression," *Advances Statist. Anal.*, vol. 95, no. 2, pp. 205-217, Jun. 2011.

[26] X. Wu, X. H. Yan, Y. H. Jo, and W. T. Liu, "Estimation of subsurface temperature anomaly in the North Atlantic using a self-organizing map neural network," *J. Atmos. Ocean. Technol.*, vol. 29, pp. 1675-1688, Nov. 2012.

[27] Y. Xu and L. Fu, "The effects of altimeter instrument noise on the estimation of the wavenumber spectrum of the sea surface height," *J. Phys. Oceanogr.*, vol. 42, no. 12, pp. 2229-2233, Dec. 2012.

<figure><img src="image_83.png" /></figure>

Pierre Tandeo was born in France in 1983. He received the M.S. degree in applied statistics from Agrocampus Ouest, Rennes, France, and the Ph.D. degree from the Oceanography from Space Laboratory, Institut Français de Recherche pour l'Exploitation de la MER, Plouzané, France, in 2010.

He spent two years as a Postdoctoral Researcher with the Atmospheric Science Research Group, University of Corrientes, Corrientes, Argentina. Since 2012, he has been a Postdoctoral Researcher with Telecom Bretagne, France. His main research interests are focused on geosciences, data assimilation, inverse problem, remote sensing data, stochastic processes, and statistical modeling.
<figure><img src="image_84.png" /></figure>

**Bertrand Chapron** was born in Paris, France, in 1962. He received the B.Eng. degree from the Institut National Polytechnique, Grenoble, France, in 1984 and the Doctorat National (Ph.D.) degree in fluid mechanics from the University of Aix-Marseille II, Marseille, France, in 1988.

He spent three years as a Postdoctoral Research Associate with the National Aeronautics and Space Administration/Goddard Space Flight Center/Wallops Flight Facility, Wallops Island, VA, USA. He has experience in applied mathematics, physical

oceanography, electromagnetic wave theory, and its application to ocean remote sensing. He is currently responsible for the Oceanography from Space Laboratory, IFREMER, Plouzané, France.

<figure><img src="image_85.png" /></figure>

**Emmanuelle Autret** was born in France in 1974. She received the M.Sc. degree in oceanography and meteorology from the University of Toulon, La Garde, France, in 2000. She is currently working toward the Ph.D. degree in signal and image processing from the Oceanography from Space Laboratory, Institut Français de Recherche pour l'Exploitation de la MER (IFREMER), Plouzané, France.

She joined IFREMER in 2002 and the Oceanography from Space Laboratory in 2005. Her research interests are ocean remote sensing, upper ocean dynamics, and data processing, particularly from high-resolution sea-surface-temperature data.

<figure><img src="image_86.png" /></figure>

**Sileye Ba** obtained the M.S. degree in mathematics, computer vision, and machine learning from the Ecole Normale Superieure de Cachan in Paris, France, in 2002.

From September 2003 to April 2009, as a Ph.D. student of the Institut de la Fondation Dalle molle d'Intelligence Artificielle Perceptive Research Institute with an affiliation to Ecole Polytechnique Federale de Lausanne (Switzerland) and a Postdoctoral Researcher of IDIAP, he worked on probabilistic methods for head pose tracking and human behavior

recognition from audio video data. From May 2009 to February 2013, he was as a Postdoctoral Researcher with the Signal and Communications Department, Telecom Bretagne, Brest, France, working on variational data assimilation methods for dynamic multimodal ocean geophysical variable modeling from multimodal satellite image sequences. Since March 2013, as a Research and Development Engineer with RN3D Innovation Lab, a start-up in Marseille (France), he has worked on computer vision and machine learning methods for near infrared image sequence analysis.

<figure><img src="image_87.png" /></figure>

**Ronan Fablet** graduated from the Ecole Nationale Supérieure de l'Aéronautique et de l'Espace (SUPAERO), Toulouse, France, in 1997. He received the Ph.D. degree in signal processing and telecommunications from the University of Rennes, Rennes, France, in 2001.

In 2002, he was an Institut National de Recherche en Informatique et Automatique Postdoctoral Fellow with Brown University, Providence, RI, USA. From 2003 to 2007, he held a full-time research position at IFREMER Brest in the field of signal and image processing applied to fisheries science. In 2008, he joined the Signal and Communications Department, Telecom Bretagne, as an Associate Professor, and has been holding a Professor position since 2012. In 2011, he was a Visiting Researcher at Institut de Recherche pour le Développement/Instituto del MAR del PErú, Peru (Peruvian Sea Research Institute). His main research interests include statistical methods for signal processing and computer vision and applications to ocean remote sensing.
4.6 Tandeo, Ailliot et Autret (2011) [SERRA]

**Contexte** Cette étude correspond au résultat principal de ma thèse, en collaboration avec l'UBO et l'IFREMER. L'article traite de la variabilité temporelle de la SST qui, à méso-échelle, explique environ 70% de la variabilité totale de cette variable océanique de surface (contre 30% pour la variabilité spatiale). Dans ce travail, nous proposons un cadre mathématique rigoureux pour la prise en compte de l'échantillonnage temporel irrégulier des observations satellitaires infrarouges, induit par la présence de nuages. Ce développement méthodologique est toujours déployé de façon opérationnelle à l'IFREMER pour la production de produits L4 de SST (V3 du produit ODYSSEA, disponible ici : https://data.marine.copernicus.eu/product/SST_ATL_SST_L4_NRT.OBSERVATIONS_010_025/description).

**Résumé** Les satellites fournissent des informations importantes pour de nombreuses variables météorologiques et océanographiques. La modélisation espace-état est couramment utilisée pour analyser de telles données avec les erreurs associées. Dans ce travail, nous proposons d'étendre le modèle espace-état linéaire Gaussien aux séries temporelles ayant un échantillonnage temporel irrégulier, comme pour celles issues de capteurs satellite. Nous discutons de l'estimation des paramètres à l'aide d'une méthode des moments et de la méthode du maximum de vraisemblance. Les résultats sur des simulations indiquent que la méthode des moments conduit à une procédure d'estimation efficace et robuste, adaptée à l'initialisation de l'algorithmes expectation-maximization, qui est combiné à une procédure d'optimisation numérique standard. La méthode proposée est ensuite validée sur des données de température de surface (SST) provenant d'un satellite infrarouge. Les résultats indiquent que la méthodologie proposée peut être utilisée pour reconstruire des séries temporelles réalistes de SST à un endroit spécifique, et qu'elle fournit des informations utiles sur la qualité des mesures satellitaires et la dynamique de la SST.
Linear Gaussian state-space model with irregular sampling:
application to sea surface temperature

Pierre Tandeo · Pierre Ailliot · Emmanuelle Autret

© Springer-Verlag 2010

**Abstract** Satellites provide important information on many meteorological and oceanographic variables. State-space models are commonly used to analyse such data sets with measurement errors. In this work, we propose to extend the usual linear and Gaussian state-space to analyse time series with irregular time sampling, such as the one obtained when keeping all the satellite observations available at some specific location. We discuss the parameter estimation using a method of moment and the method of maximum likelihood. Simulation results indicate that the method of moment leads to a computationally efficient and numerically robust estimation procedure suitable for initializing the Expectation-Maximisation algorithm, which is combined with a standard numerical optimization procedure to maximize the likelihood function. The model is validated on sea surface temperature (SST) data from a particular satellite. The results indicate that the proposed methodology can be used to reconstruct realistic SST time series at a specific location and also give useful information on the quality of satellite measurement and the dynamics of the SST.

**Keywords** State-space model · Irregular sampling · Ornstein–Uhlenbeck process · EM algorithm · Sea surface temperature

# 1 Introduction

Sea surface temperature (SST) is an important oceanographic variable for many applications (see e.g. [7] and references therein). Several satellites and buoy networks provide continuous observations of this variable, leading to a huge amount of data. Statistical methods are then needed to combine all this information and provide realistic SST analysis at any date and any location in the ocean.

State-space models provide a flexible methodology for analysing such complex environmental data sets, and they have already been used in a wide range of problems (see e.g. [13]) including meteorological and oceanographic applications (see e.g. [1, 11, 17, 26]). The basic idea of these models consists in introducing the “true” value of the physical variable of interest as a hidden variable (the “state”). Then, stochastic models are used both to describe the dynamics of the state and to relate the observations to the state. When linear Gaussian models are used, we get the so-called linear Gaussian state-space model which has been extensively studied in the literature (see e.g. [8] and references therein). Note that [14] proposed unified notations for state-space models and data assimilation in oceanography and meteorology which are partially adopted here.

In this work we analyse satellite SST data at a single location, where buoy data is available for comparison, and we consider the time series obtained by keeping all the satellite data available nearby this location. It leads to a time series with irregular time-step, with generally several data each day but also sometimes gaps of several days with

<footnote>P. Tandeo (✉) · E. Autret
Laboratoire d'Océanographie Spatiale, IFREMER, Plouzané,
France</footnote>

<footnote>e-mail: pierre.tandeo@ifremer.fr</footnote>

<footnote>E. Autret
e-mail: emmanuelle.autret@ifremer.fr</footnote>

<footnote>P. Ailliot
Laboratoire de Mathématiques, UMR 6205, Université
Européenne de Bretagne, Brest, France</footnote>

<footnote>e-mail: pierre.ailliot@univ-brest.fr</footnote>
no data. We adopt a continuous-time state-space model to analyse this time series in which the state is supposed to be an Ornstein–Uhlenbeck process. It leads to a simple generalization of the usual linear Gaussian state-space model with regular time-step.

The most usual method for estimating the parameter in models with latent variable consists in computing the maximum likelihood estimates using the Expectation-Maximisation (EM) algorithm. In this work, we propose to improve the numerical efficiency of the EM algorithm by combining it with a method of moment and a standard numerical optimization procedure. The method of moment is used to provide realistic starting values to the EM algorithm with the extra benefit of providing graphical tools which permit to asses the realism of the model. The standard numerical optimization procedure is used to accelerate the convergence of the EM algorithm near the maxima and provide estimates of the observed information matrix and thus important information on the variance of the estimates.

The paper is organised as follows. The SST data and the model are introduced in Sect. 2. Then, the parameter estimation is discussed in Sect. 3: after describing the practical implementation of the various methods, we asses the efficiency of the whole procedure through simulations. In Sect. 4, we discuss the results obtained on the data with the proposed methodology. Conclusions are drawn in Sect. 5.

## 2 Data and model

Several instruments on-board satellites provide measurements of SST over the entire surface of the ocean with different spatial and temporal resolutions. In this work, we focus on the data provided by the infrared Advanced Very High Resolution Radiometer instrument on-board the METOP satellite (see [18] for more details). This satellite covers the global ocean with a spatial resolution of 0.05° and provides two SST observations per day at the most in

optimal conditions. In this paper, we first consider the data available at a given location, with geographical coordinates (0°N, 23°W), in the tropical region of the Atlantic Ocean. More precisely, we consider 2 years of data, from 11-Jul-2007 to 18-Jun-2009, which are representative of the variability of the SST conditions at this location. Hereafter, ($t_1, ..., t_n$) denotes the times at which the METOP satellite data are available, with $n = 1087$ the total number of observations. Since satellite observations may be contaminated by atmospheric conditions (e.g. cloud coverage), some data are missing and the time difference $\Delta_i = t_i - t_{i-1}$ between two consecutive observations may vary from a half day to a few days (see Fig. 1).

The resulting time-series is clearly non-stationary (see Fig. 2) with in particular important seasonal components. The non-stationary components have complex features and we could not find any appropriate parametric model to describe them. We have thus decided to use the SST analysis produced by the National Climatic Data Center (NCDC) (daily “OIV2 analysis” with 0.25° spatial resolution) to remove these components. These analysis are derived from different satellite sources independent of METOP data (see [22]) and we assume that they provide a good estimate of the low-variations of the SST conditions. Both data sources METOP and OIV2 are available at the URL http://www.hrdds.net.

Then we consider the time series $y_{t1}^{t_n} = (y_{t_1}, ..., y_{t_n})$ obtained by removing the OIV2 analysis from the METOP data (see Fig. 3). We assume that this new time-series, referred as the SST anomaly hereafter, is a discrete-time realization of a continuous-time stationary process {$Y_t$}. Modelling the time series {$y_{t1}^n$} may provide important information on the small scale variability of SST and also on the quality of METOP measurements and OIV2 analysis as it will be shown in Sect. 4 and finally lead to a better assimilation of these data into numerical models.

The model that we consider for {$Y_t$} is introduced below. First, we assume that the observed SST anomaly at

<figure><img src="image_88.png" /><figcaption>**Fig. 1** Histogram of the time lags $\Delta_t$ in days (**left**) and of the SST anomalies {$y_{t1}^n$} in °C (**right**)</figcaption></figure>
<figure><img src="image_89.png" /><figcaption>**Fig. 3** SST anomalies (in °C) obtained by removing the OIV2 analysis from METOP data</figcaption></figure>

<figure><img src="image_90.png" /></figure>

time $t$, $Y_t$, is related to the “true” SST anomaly at time $t$, denoted $X_t$, by the measurement equation below:

$$Y_t = HX_t + \sqrt{R}\epsilon_t \quad (1)$$

where {$\epsilon_t$} is a Gaussian white noise sequence with zero mean and unit variance. In practice $R$ represents the variance of the observation error and $H$ allows a transformation between the state and the observations. For the particular METOP measurements considered in this paper (we keep only the best quality data), the standard deviation of the observation error has been estimated globally to 0.5°C, but it is known that it may vary according to the retrieval algorithm (day-time and night-time), the region and the season (see [18] for more details). The observation equation (1) could be modified to take into account these fluctuations in the accuracy of the data. In the same way, we could include the various covariates which alter the quality of the satellite measurements (see [25]) or assume that the parameters $H$ and $R$ depend on the satellite if the observed time series was obtained by mixing data from different satellites.

Then we assume that the latent process {$X_t$} is a simple Ornstein–Uhlenbeck process, that is a stationary solution of the following stochastic differential equation:

$$dX_t = -\lambda X_t dt + \tau dW_t \quad (2)$$

where {$W_t$} denotes a standard Brownian motion. A physical justification of using this model to describe the local dynamics of the SST, when neglecting horizontal transport and heat exchange, is given in [10]: $\lambda > 0$ is the time correlation (in day) or feedback parameter which represents the slowly evolving transfer of heat and $\tau > 0$ the variability coming from weather fluctuations (see also [19, 21]).

Hereafter, we denote $\sigma^2 = \text{Var}(X_t) = \frac{\tau^2}{2\lambda}$ the variance of the stationary distribution. {$X_t$} is a Markov process which satisfies, for $i \in \{2, \dots, n\}$,

$$X_{t_i} = M_{\Delta_i} X_{t_{i-1}} + \sqrt{Q_{\Delta_i}} \eta_{t_i} \quad (3)$$

with $M_{\Delta_i} = \exp(-\lambda \Delta_i)$, $Q_{\Delta_i} = \sigma^2(1 - M_{\Delta_i}^2)$ and {$\eta_{t_i}$} $i \in \{2, \dots, n\}$ a Gaussian white noise sequence with zero mean and unit variance independent of {$\epsilon_{t_i}$} $i \in \{1, \dots, n\}$. In the particular case when the temporal sampling is regular, i.e. when $\Delta_1 = \dots = \Delta_n$, we retrieve a standard AR(1) process and the usual linear Gaussian state-space model. Here again, more complicated models could be considered, with for example non-linear dynamics, but this would complicate the statistical inference methods discussed in the next section.

Finally, the various conditional independence assumptions which imply the particular Markovian structure of the
<table><caption>Fig. 4 Directed acyclic graph for the linear Gaussian state-space model with irregular time step</caption><tr><td>State</td><td>x<sub>t<sub>i-1</sub></sub><br>↓</td><td>→</td><td>x<sub>t<sub>i</sub></sub><br>↓</td><td>→</td><td>x<sub>t<sub>i+1</sub></sub><br>↓</td><td>...</td></tr><tr><td>Observation</td><td>y<sub>t<sub>i-1</sub></sub></td><td></td><td>y<sub>t<sub>i</sub></sub></td><td></td><td>y<sub>t<sub>i+1</sub></sub></td><td></td></tr><tr><td>Time</td><td>t<sub>i-1</sub><br>...</td><td>...</td><td>t<sub>i</sub><br>Δ<sub>i</sub></td><td>...</td><td>t<sub>i+1</sub><br>Δ<sub>i+1</sub></td><td>...</td></tr></table>

state-space model, when observed at discrete time $t_1, \dots, t_n$, are summarized on the directed acyclic graph shown on Fig. 4.

## 3 Parameter estimation

The estimation of the unknown parameters in Gaussian linear state-space models observed at regular time step has been addressed by many authors and the most usual method consists probably in computing the maximum likelihood (ML) estimates using the EM algorithm (see e.g. [8]).

However, before computing the ML estimates, it is important to check the identifiability of the parameters. For the particular model under consideration, it is possible to show that the observations follow a multivariate Gaussian distribution with an explicit covariance function. Using this result, we can give conditions on the parameters which ensure identifiability and also propose a first method based on the moments to estimate the parameters. The corresponding estimates will be denoted MOM estimates hereafter. This is discussed in Sect. 3.1 Then, in Sect. 3.2, we detail the practical implementation of the EM algorithm for the Gaussian linear state-space model with irregular time-step. We discuss how it can be combined with the method of moment and a more standard numerical optimization procedure proposed in [16] to get a computationally efficient and numerically robust estimation procedure. Finally, this is illustrated in Sect. 3.3 through simulations.

### 3.1 Covariance function

With the various assumption made in the previous section, {$Y_t$} is a stationary Gaussian process with zeros mean and covariance function

$$ \mathrm{Cov}(Y_t, Y_{t'}) = H^2 \sigma^2 \exp(-\lambda |t - t'|) + R \mathbf{1}_{\{0\}}(t - t') \quad (4) $$

We deduce that the distribution of the observed sequence ($y_{t_1}, \dots, y_{t_n}$) is a multivariate Gaussian distribution with zeros mean and covariance matrix which can be expressed from the unknown parameter $H$, $R$, $\sigma^2$ and $\lambda$. According to (4), this covariance matrix depends on the parameters $H$ and $\sigma^2$ only through the product $H^2\sigma^2$ and thus we need

to add a constraint in order to ensure identifiability of the parameters. Hereafter, we fix $H = 1$ and denote $\theta = (\lambda, \sigma^2, R) \in (0, +\infty)^3$ the unknown parameters.

The covariance function (4) corresponds to a classical model in spatial statistics since we retrieve an exponential model with nugget $R$, sill $\sigma^2 + R$ and range $1/\lambda$. Usual methods in geostatistics permit to compute an empirical estimate of the variogram from the data (see e.g. [3, p. 69]). The variogram is directly related to the covariance function for second order stationary processes and the empirical variogram can be used to check the realism of the parametric model 4 and also fit it using the weighted least square method. Here the weights depend on the number of pairs of time points which are available to estimate the empirical variogram as discussed in [3, p. 96]. The corresponding estimates will be denoted MOM estimates hereafter.

### 3.2 Maximum likelihood estimation

Alternatively, the parameters can be estimated by computing the ML estimates. According to the conditional independence assumptions shown on Fig. 4, the complete log-likelihood, based on both the latent and observed sequences, is given by

$$ \log(p(x_{t_1}^{t_n}, y_{t_1}^{t_n}; \theta)) = \log(p(x_{t_1})) + \sum_{i=2}^{n} \log(p(x_{t_i} | x_{t_{i-1}}; \theta)) \\ + \sum_{i=1}^{n} \log(p(y_{t_i} | x_{t_i}; \theta)) $$

where the conditional distributions $p(x_{t_i}|x_{t_{i-1}}; \theta)$ and $p(y_{t_i}|x_{t_i}; \theta)$ are Gaussian distributions which characteristics are given respectively by (3) and (1). Hereafter we will assume that the initial distribution $p(x_{t_1})$ is a Gaussian distribution with known mean $x^{(b)}$ and variance $B$ and in practice these values will be estimated using historical data. Thus, apart from a constant, we obtain

$$ \log(p(x_{t_1}^{t_n}, y_{t_1}^{t_n}; \theta)) = -(n-1) \log(\sigma) - \frac{1}{2} \sum_{i=2}^{n} \log(1 - \exp(-2\lambda\Delta_i)) \\ - \frac{1}{2\sigma^2} \sum_{i=2}^{n} \frac{(x_{t_i} - \exp(-\lambda\Delta_i)x_{t_{i-1}})^2}{(1 - \exp(-2\lambda\Delta_i))} \\ - \frac{n}{2} \log(R) - \frac{1}{2R} \sum_{i=1}^{n} (y_{t_i} - x_{t_i})^2 \quad (5) $$

The ML estimates $\hat{\theta}$ is the value of $\theta$ that maximises the (incomplete) likelihood of the observations $y_{t_i}^{t_n}$ formed by integrating the complete likelihood (5) over the missing variables.
In this paper, the EM algorithm due to [4] is used to compute $\hat{\theta}$. This recursive algorithm computes successive approximations $\hat{\theta}_k = (\lambda_k, \sigma_k^2, R_k)$ of $\hat{\theta}$ by cycling through the following steps.

**E-step:** Compute $U(\theta|\hat{\theta}_k) = E(\log(p(X_{t_1}^{t_n}, y_{t_1}^{t_n}; \theta))|y_{t_1}^{t_n}; \hat{\theta}_k)$ as a function of $\theta$.

**M-step:** Determine the updated parameter estimate $\hat{\theta}_{k+1} = \arg\max_{\theta} U(\theta|\hat{\theta}_k)$.

Under certain general conditions it can be shown that the sequence of estimates $\hat{\theta}_n$ yields monotonically increasing values of the incomplete likelihood, and converges to a maximum of this function (see [15]). Thus the EM algorithm provides an alternative method of maximising the incomplete log-likelihood which is commonly used in models with hidden or latent variables such as the model proposed here. The EM algorithm directly utilises the hidden structure and, as a consequence, is often more robust in practice to the choice of starting values than usual numerical optimization methods. Its computational efficiency is enhanced if the E and M steps are readily evaluated. Various authors have discussed the practical implementation of these steps for linear Gaussian state-space models with regular time sampling ([2, 5, 6, 24], pp. 384–388). Hereafter, we discuss the extension to the case with irregular sampling.

**E step** To determine $U(\theta|\hat{\theta}_k)$ as a function of $\theta$ we need to compute the following smoothing probabilities, for $i=1, \dots, n$:

$$ 
\begin{aligned}
x_{t_i}^{(s)} &= E(X_{t_i}|y_{t_i}^{t_n}; \hat{\theta}_k), & x_{t_i,t_i}^{(s)} &= E(X_{t_i}^2|y_{t_i}^{t_n}; \hat{\theta}_k), \\
x_{t_{i-1},t_i}^{(s)} &= E(X_{t_{i-1}}|X_{t_i}y_{t_i}^{t_n}; \hat{\theta}_k)
\end{aligned} 
\quad (6) $$

These quantities can be computed using the Kalman recursions described hereafter. This is a particular case of the general Kalman recursions given for example in [24] and [2] pp. 127–147.

**• Kalman filter** Let us denote

$$ x_{t_i}^{(f)} = E(X_{t_i}|y_{t_i}^{t_{i-1}}; \hat{\theta}_k), \quad P_{t_i}^{(f)} = \text{Var}(X_{t_i}|y_{t_i}^{t_{i-1}}; \hat{\theta}_k) $$

the mean and the variance of the forecast probabilities and

$$ x_{t_i}^{(a)} = E(X_{t_i}|y_{t_i}^{t_i}; \hat{\theta}_k), \quad P_{t_i}^{(a)} = \text{Var}(X_{t_i}|y_{t_i}^{t_i}; \hat{\theta}_k) $$

the mean and the variance of the filtering probabilities. These quantities can be computed using the recursion below.

**Initialization:** Compute the Kalman filter gain $K_{t_i} = \frac{B}{B+R}$ and

$$ x_{t_1}^{(a)} = x^{(b)} + K_{t_1}(y_{t_1} - x^{(b)}), \quad P_{t_1}^{(a)} = (1 - K_{t_1})B $$

where the parameters $x^{(b)} = E[X_{t_1}]$ and $B = \text{Var}(X_{t_1})$ of the initial distribution are supposed to be known.

**Recursion:** for $i=2, \dots, n$

**- Time update:**

$$ x_{t_i}^{(f)} = M_{\Delta_i} x_{t_{i-1}}^{(a)}, \quad P_{t_i}^{(f)} = M_{\Delta_i}^2 P_{t_{i-1}}^{(a)} + Q_{\Delta_i} $$

**- Observation update:** compute the Kalman filter gain

$$ K_{t_i} = \frac{P_{t_i}^{(f)}}{P_{t_i}^{(f)} + R} \quad \text{and} $$

$$ x_{t_i}^{(a)} = x_{t_i}^{(f)} + K_{t_i} (y_{t_i} - x_{t_i}^{(f)}), \quad P_{t_i}^{(a)} = (1 - K_{t_i}) P_{t_i}^{(f)} $$

**• Kalman smoother** Let us denote

$$ P_{t_i}^{(s)} = \text{Var}(X_{t_i}|y_{t_i}^{t_n}; \hat{\theta}_k) $$

the variance of the smoothing probabilities at time $t_i$. These quantities and the conditional expectation $x_{t_i}^{(s)}$ define in (6) can be computed using the backward recursions below.

**Initialization:**

$$ x_{t_n}^{(s)} = x_{t_n}^{(a)}, \quad P_{t_n}^{(s)} = P_{t_n}^{(a)} $$

**Recursion:** for $i=n-1, \dots, 1$ compute the Kalman smoother gain $K_{t_i}^{(s)} = \frac{P_{t_i}^{(a)}M}{P_{t_{i+1}}^{(f)}}$ and

$$ 
\begin{aligned}
x_{t_i}^{(s)} &= x_{t_i}^{(a)} + K_{t_i}^{(s)} (x_{t_{i+1}}^{(s)} - x_{t_{i+1}}^{(f)}), \\
P_{t_i}^{(s)} &= P_{t_i}^{(a)} + (K_{t_i}^{(s)})^2 (P_{t_{i+1}}^{(s)} - P_{t_{i+1}}^{(f)}),
\end{aligned} 
$$

Finally $U(\theta|\hat{\theta}_k)$ can be computed from the quantities computed with the Kalman smoother above and the relations

$$ 
\begin{aligned}
x_{t_i,t_i}^{(s)} &= P_{t_i}^{(s)} + (x_{t_i}^{(s)})^2, \\
x_{t_{i-1},t_i}^{(s)} &= \text{Cov}(X_{t_{i-1}}, X_{t_i}|y_{t_i}^{t_n}; \hat{\theta}_k) + x_{t_{i-1}}^{(s)}x_{t_i}^{(s)}
\end{aligned} 
$$

where

$$ 
\begin{aligned}
\text{Cov}(X_{t_{i-1}}, X_{t_i}|y_{t_i}^{t_n}; \hat{\theta}_k) &= (1 - K_{t_i})MP_{t_{i-1}}^{(a)} \\
&\quad + \frac{P_{t_i}^{(s)} - P_{t_i}^{(a)}}{P_{t_{i-1}}^{(a)}}(1 - K_{t_i})MP_{t_i}^{(a)}
\end{aligned} 
$$

**M step**

The function $U(\theta|\hat{\theta}_k)$ can be decomposed as

$$ U(\theta|\hat{\theta}_k) = U_X(\lambda, \sigma^2|\hat{\theta}_k) + U_{Y|X}(R|\hat{\theta}_k) $$

where

$$ 
\begin{aligned}
U_X(\lambda, \sigma^2|\hat{\theta}_k) \\
&= -(n-1)\log(\sigma) - \frac{1}{2}\sum_{i=2}^n \log(1 - \exp(-2\Delta_i\lambda)) \\
&\quad - \frac{1}{2\sigma^2}\sum_{i=2}^n \frac{x_{t_i,t_i}^{(s)} - 2\exp(-\Delta_i\lambda)x_{t_{i-1},t_i}^{(s)} + \exp(-2\Delta_i\lambda)x_{t_{i-1},t_{i-1}}^{(s)}}{1 - \exp(-2\Delta_i\lambda)}
\end{aligned} 
$$
and

$$U_{Y|X}(R|\hat{\theta}_k) = -\frac{n}{2}\log(R) - \frac{1}{2R}\sum_{i=1}^{n}\left\{y_{t_i}^2 - 2y_{t_i}x_{t_i}^{(s)} + x_{t_i,t_i}^{(s)}\right\}$$

The second term $U_{Y|X}$ is similar to the case with regular sampling and the maximum is obtained for $R = R_{k+1}$ with

$$R_{k+1} = \frac{1}{n}\sum_{i=1}^{n}\left\{y_{t_i}^2 - 2y_{t_i}x_{t_i}^{(s)} + x_{t_i,t_i}^{(s)}\right\}$$

The first term $U_X$ is specific to the case with irregular sampling and numerical optimisation procedures have been used to compute $(\lambda_{k+1}, \sigma_{k+1}^2)$ since we could not derive analytic expressions these quantities. Here the relation

$$\sigma_{k+1}^2 = \frac{1}{n-1}\sum_{i=2}^{n} \frac{x_{t_i,t_i}^{(s)} - 2\exp(-\Delta_i\lambda_{k+1})x_{t_{i-1},t_i}^{(s)} + \exp(-2\Delta_i\lambda_{k+1})x_{t_{i-1},t_{i-1}}^{(s)}}{1 - \exp(-2\Delta_i\lambda_{k+1})}$$

has been used to transform the initial two-dimensional optimization problem into a simple one-dimensional optimisation problem and reduce computational time.

The EM algorithm has several well known limitations. First it may converge to a non-interesting local maximum of the likelihood function depending on the starting value $\hat{\theta}_0$, and thus it is important to provide realistic initial parameter values. Here we have used the estimates obtained using the method of moment described in Sect. 3.1 Indeed the various tests that we have done indicate that this method leads to robust estimates and generally provide a good starting value to the EM algorithm with low numerical cost (see Sect. 3.3). This is particularly useful to avoid numerical problem when fitting the model to a large number of data sets for regional studies such as the one performed in Sect. 4.4

Another limitation of the EM algorithm is its slow convergence near the maxima where using a standard optimization algorithm is generally far more efficient, at least when it is possible to compute the incomplete likelihood function quickly. For the model under consideration, the incomplete likelihood function is a sub-product of the Kalman filter since we have

$$p(y_{t_i}^{d_i}; \theta) = \prod_{i=2}^{n} p(y_{t_i} | y_{t_{i-1}}^{d_{i-1}})$$

where the conditional distribution $p(y_{t_i}|y_{t_{i-1}}^{d_{i-1}})$ is a Gaussian distribution with mean $E(X_{t_i}|y_{t_{i-1}}^{d_{i-1}})$ and variance $Var(X_{t_i}|y_{t_{i-1}}^{d_{i-1}}) + R$ and these quantities are computed recursively in the Kalman filter (see Sect. 3.2). Eventually, the gradient of the log-likelihood function could also be computed to accelerate the convergence of the numerical

optimization procedure. In this work, we did not provide the gradient to the Matlab function used for the numerical optimization but we did not encounter any numerical problem and the computational efficiency was good enough.

Another advantage of switching from the EM algorithm to a quasi-Newton algorithm close to the maxima is that quasi-Newton algorithms provide an approximation of the Hessian of the log-likelihood function, and thus useful information on the variance of the ML estimates (see Sect. 3.3)

### 3.3 Simulations

In this section, the relative performances of ML and MOM estimates are assessed through simulations. More precisely, for various values of $n \in \{200, 300, \dots, 2000\}$, we have simulated $N = 1,000$ sequences of length $n$ using the scheme described below:

1. Simulate the time lags $(\Delta_i)_{i \in \{2, \dots, n\}}$ as an i.i.d. sample from the empirical distribution of the time lags for satellite data (see Fig. 1).

2. Simulate the initial state $x_{t_1}$ as a Gaussian variable with mean $x^{(b)}$ and variance $B$ and then recursively $(x_{t_i})_{i \in \{2, \dots, n\}}$ according to (3).

3. Simulate the observed process $(y_{t_i})_{i \in \{1, \dots, n\}}$ using (1).

The following parameters values have been chosen for the numerical experiment: $\lambda = 0.5$, $B = \sigma^2 = 0.05$, $R = 0.5$ and $x^{(b)} = 0$. It corresponds to realistic values for the application discussed in the next section.

Then, for each simulated sequence the ML and MOM estimates have been computed. In practice, ML and MOM estimates have been computed using a quasi-Newton algorithm with the true values of the parameters as initial value. Although such initialization is not possible for practical applications, it permits to avoid convergence to non interesting local maxima of the likelihood function and a fair comparison of the two estimates. Figure 5 shows the empirical estimate of the bias and variance of the estimates computed from these simulations. As expected, the ML estimates generally outperform the MOM estimates in terms of both bias and variance. However, the MOM estimates give satisfactory results for the different values of $n$ and have the advantage of being computed with low computational costs and less sensitive to the choice of realistic starting values than the EM algorithm. For comparison purpose, the variances computed from the inverse of the observed information matrix are also shown on Fig. 5. The agreement with the empirical variances of the ML estimates is generally good, especially for large sample size as expected from the general asymptotic theory for the ML estimates.
<figure><img src="image_91.png" /><figcaption>**Fig. 5** Plot of the simulated bias (**left**) and variances (**right**) for the MOM (**dashed-dotted line**) and ML estimates (**full line**) for different length sequences **n** (x axis). Estimate of $\lambda$ (**top panel**), of $\sigma^2$ (middle) and **R** (**bottom**). The **dotted lines** on the **right panel** is the variance computed from the information matrix (empirical mean over the different simulations). The simulated results are based on $N = 1,000$ replications</figcaption></figure>

# 4 Application to SST data

In this section the model is first fitted and validated on the SST data introduced in Sect. 2. The original time series has been divided into two consecutive parts: the first one ($y_{t1}, \dots, y_{t_{n_1}}$) for estimating the parameters and second one ($y_{t_{n_1+1}}, \dots, y_{t_n}$) for validating the model. In practice, we used $n_1 = 725$ observations to fit the model, a reasonable amount of data according to the simulation results given in Sect. 3.3 It corresponds to a proportion of about two-thirds of the data (more than 1 year).

In Sect. 4.1, we first discuss the results obtained when fitting the model on the training data set. Then the model is validated using cross-validation on the validation data set in Sect. 4.2 and by comparison to buoy data in Sect. 4.3. Finally, in Sect. 4.4, the methodology is applied to data at many locations on a regular grid covering the Atlantic ocean and the spatial behaviour of the parameter estimates is discussed.

## 4.1 Parameter estimation

The parametric covariance model (4) has been fitted to the empirical estimate of the autocovariance function of the

SST anomaly using weighted least square method leading to the MOM estimates (see Sect. 3.1). The corresponding variograms are shown in Fig. 6. The overall agreement is good, except maybe a 5 days component which is visible on the empirical variogram function (see [12] for a discussion on the existence of peak frequencies in SST time series). This indicates that the assumptions made on the shape of the covariance function is realistic, at least when focussing to time lags up to 40 days. Let us remark that according to Fig. 1, it seems also reasonable to assume that the marginal distribution is approximately Gaussian except maybe the lower tail of the distribution.

Starting from the MOM estimates obtained by fitting the covariance function, we have run the EM algorithm. The first iterations are efficient and the likelihood function increases rapidly (see Fig. 7) but after some iterations the convergence becomes rather slow, and switching to a standard numerical optimisation procedure permits to save computational time. According to Table 1, the ML estimate of $\lambda$ is significantly lower than the MOM estimate and the ML estimates of $\sigma^2$ and $R^2$ are higher than the corresponding MOM estimates, although the differences for $\sigma^2$ and $R^2$ do not seem to be statistically significant if we compare the differences in the parameter values to the
<figure><img src="image_92.png" /><figcaption>**Fig. 7** Evolution of the parameters values during the 100 iterations (x axis) of the EM algorithm: $\hat{\lambda}$ (top-left), $\hat{\sigma}^2$ (top-right), $\hat{R}$ (bottom-left). The bottom-right panel shows the increase of the log-likelihood function</figcaption></figure>

<figure><img src="image_93.png" /><figcaption>**Table 1** Parameter value after the different steps of the fitting procedure: method of moment (first column), 100 iterations of the EM algorithm (second column) and numerical optimization of the likelihood function with a quasi-Newton algorithm (third column)</figcaption></figure>

<table><thead><tr><th rowspan="2"></th><th rowspan="2">Method of moments</th><th colspan="2">Maximum likelihood</th><th rowspan="2">Standard deviation</th></tr><tr><th>EM algorithm</th><th>Quasi-Newton</th></tr></thead><tbody><tr><td>$\hat{\lambda}$ (day<sup>-1</sup>)</td><td>0.145</td><td>0.057</td><td>0.056</td><td>0.019</td></tr><tr><td>$\hat{\sigma}^2$</td><td>0.184</td><td>0.329</td><td>0.330</td><td>0.094</td></tr><tr><td>$\hat{R}$</td><td>0.129</td><td>0.141</td><td>0.141</td><td>0.010</td></tr><tr><td>Log-likelihood</td><td>-463.15</td><td>-450.69</td><td>-450.68</td><td></td></tr></tbody></table>

The last column gives an estimate of the standard deviation of the ML estimates computed from the information matrix. Results obtained on the training data set

standard deviations given in Table 1. ML estimates identify a second-order structure with a higher sill, which better coincides with the empirical variance of the time series (about 0.47), and also a higher range. Despite these differences in the parameters values, the agreement between the covariance functions is good for time lags less than 10 days (see Fig. 6) and thus we may expect that we would get similar results if using the model with the MOM instead of the ML estimates for estimating the true SST in Sects. 4.2 and 4.3.

The final parameter values are in good agreement with our knowledge of the physical process under consideration. In particular, according to [18], the standard deviation of the measurement error of the METOP data considered in this paper may vary between 0.33 and 0.51°C depending on the conditions. This range matches with the 95% confidence interval for the standard deviation of the observation error (we get approximately the interval between 0.35 and 0.40°C). Then, the low value of $\hat{\lambda}$ imply an important temporal persistence of the SST conditions and is coherent with the climatology of the place of interest were SST anomaly is known to have a strong temporal correlation. Finally, comparing the variance of the innovation of the dynamics for a time lag of 1 day ($\hat{Q}_1 = 0.04$) with the one of observation error $\hat{R}$ indicates that more weights will generally be given to the previous analysis than to the current observation in the Kalman recursions.

## 4.2 Cross-validation

In this section, we validate the model using cross-validation on the validation data set. For each $i \in \{n_1 + 1, \dots, n\}$, the observation at time $t_i$ is removed and the Kalman recursions are used to compute
$$ 
\begin{aligned}
x_{t_i|t}^{(s)}(\hat{\theta}) &= E(X_{t_i} | y_{t_{n+1}}, y_{t_{i+1}}^{t_n}; \hat{\theta}), P_{t_i|t}^{(s)}(\hat{\theta}) \\
&= \mathrm{Var}(X_{t_i} | y_{t_{n+1}}, y_{t_{i+1}}^{t_n}; \hat{\theta})
\end{aligned}
$$

If the various assumptions made in Sect. 2 are valid, then the conditional distribution of $Y_{t_i}$ given the past observations $y_{t_{n+1}}^{t_n}$ and the future observation $y_{t_{i+1}}^{t_n}$ should be approximately Gaussian with mean $x_{t_i|t}^{(s)}(\hat{\theta})$ and variance $P_{t_i|t}^{(s)}(\hat{\theta}) + \hat{R}$. An histogram of the standardized residuals

$$ \frac{y_{t_i} - x_{t_i|t}^{(s)}(\hat{\theta})}{\sqrt{P_{t_i|t}^{(s)}(\hat{\theta}) + \hat{R}}} $$

is shown on Fig. 8 together with the probability density function of the standard normal distribution and a normal probability plot (formal goodness of fit test are hard to implement since the residuals are not independent). The fit is generally good except again for the lower part of the distribution and this indicates that there are too many low residuals. According to Fig. 9, it corresponds to breaks in the observed time series at date when the SST anomaly suddenly drops. It is known that various factors (aerosol optical depth, wind speed or proximity to clouds for example) may perturb the quality of the data and a careful examination of these factors at the dates when the SST

drops has been done. We could not identify anything special at these dates and thus we believe that the drops are due to non-linearities in the dynamics of the true SST anomaly. It indicates that using a non-linear model instead of (2) may be more appropriate. Let us remark that the standardized residuals may also provide useful information on outliers.

## 4.3 Comparison with buoy data

Using the model proposed in this work and the Kalman smoother on SST anomaly derived from satellite data, we can estimate the "true" SST anomaly at any time and thus emulate a virtual buoy. In order to check the realism of such virtual buoy, we have been compared the result with SST buoy measurements available at high temporal resolution (10 min) from the Pilot Research Moored Array in the Tropical Atlantic (PIRATA, see [23]) at the same location (0°N, 23°W). According to Fig. 10, the virtual buoy obtained by smoothing satellite data has some similarities with buoy data, but there are also important differences (only 63% of buoy measurements are contained in the 95% fluctuation intervals for the smoothing probabilities). However, Table 2 indicates that using the model proposed in this paper permits to improve the quality of the

<figure><img src="image_94.png" /><figcaption>**Fig. 9** Raw (dotted line) and interpolated (full line) satellite SST anomalies (in °C) together with a 95% fluctuation interval for the smoothing probabilities (grey). Results obtained by cross-validation on the validation data set</figcaption></figure>

<figure><img src="image_95.png" /></figure>
<figure><img src="image_96.png" /><figcaption>**Table 2** Difference between satellite (raw and smoothed) and buoy SST (bias, standard deviation and root mean square error) computed on the validation data set</figcaption></figure>

<table><thead><tr><td></td><td>Bias</td><td>Standard deviation</td><td>RMSE</td></tr></thead><tbody><tr><td>Raw satellite data</td><td>-0.22</td><td>0.47</td><td>0.52</td></tr><tr><td>Smoothed satellite data</td><td>-0.22</td><td>0.31</td><td>0.38</td></tr></tbody></table>

original satellite data and decrease the standard deviation of the error but can not correct the negative bias present in the original satellite data (underestimation of the SST measured at the buoy).

Since the results given in the previous sections indicate that the state-space model proposed in this paper is realistic for satellite data, we may conclude that the significant differences between the buoy and the virtual buoy are due to differences in the satellite and buoy data. A first reason may be the well known depth-to-skin bias discussed in [18]: METOP satellite measures the skin SST (the temperature of the sea in the first µm) whereas the buoy measures the temperature at a depth of about 1 meter and the temperature gradient evolves strongly in this surface layer. A second possible reason is the difference in the scale of the measurements: buoy data are local measurements and are able to identify small scale variation whereas METOP data describes larger scale variations since they retrieve the mean SST over a 5 × 5 km² surface.

Finally, these results highlight the difficulties of building a realistic SST time series from satellite data. Possible improvements are discussed in the conclusion.

## 4.4 Generalization to the Atlantic ocean

The methodology introduced above for the point with geographical coordinates (0°N, 23°W) has been applied to locations on a regular grid with 1° resolution in both latitude and longitude covering the Atlantic Ocean. The state-space model is fitted at each point on the time series of SST

anomalies obtained by removing OIV2 analysis from METOP data. The length *n* of the time series depends on the location of interest and varies from 100 to 900 (see Fig. 11). According to the simulation results given in Sect. 3.3, this may lead to estimates with high variance at locations with poor satellite coverage.

The spatial behaviour of the parameter estimates shown on Fig. 11 gives important information on the small-scale variability of SST and also on the quality of METOP data and OIV2 analysis. First, the feedback parameter $\lambda$ (expressed in day$^{-1}$) informs us about the heat transfer at the surface of the ocean. In order to facilitate the interpretation, we have chosen to represent the spatial evolution of $M_1 = \exp(-\lambda)$ which corresponds to the autoregression coefficient for a time lag of 1 day between two observations. The estimate of $M_1$ mainly depends on the latitude with longer range temporal dependence in the inter-tropical convergence zone (ITCZ) than in the mid-latitudes. Then, the variance of the stationary distribution of the state $\sigma^2$ informs us about the variability of the SST anomaly. According to Fig. 11, the areas with high variability correspond to places, like the Falkland area off the Brazilian coast and the Gulf Stream off the Canadian coast, with strong sea-surface currents and wind conditions. Moreover, the more important upwelling systems of the Atlantic ocean can also be identified, e.g. the Canary and Benguela regions which are areas with strong winds yielding to a mixing of the ocean layer. In the rest of the Atlantic ocean, the variance is about 0.1. Finally, the value of the parameter *R* is the variance of the measurement errors of the METOP sensor. Estimate of this variance were provided in a previous study [18] by comparing METOP observations to data from drifting buoys. Unfortunately, the number of buoys is limited and covers a small part of ocean. The approach presented in this paper, based only on remotely sensed data, presents a global view of the spatial distribution of *R*. According to Fig. 11, the principal sources of contamination of METOP infra-red sensor seem to be the aerosol of the Saharan dust (see [9]) and the wildfire off the Angola coast.
<figure><img src="image_97.png" /><figcaption>**Fig. 11** *Top-left* number of METOP data at each grid point. *Top-right* spatial evolution of the estimate of the 1 day autocorrelation coefficient</figcaption></figure>

$\hat{M}_1 = \exp(-\hat{\lambda})$. *Bottom-left* spatial evolution of the estimate of the variance of the stationary distribution $\hat{\sigma}^2$. *Bottom-right* spatial evolution of the estimate of the variance of the measurement error *R*

# 5 Conclusion and perspectives

In this paper, we propose an extension of the usual linear and Gaussian state space model to analyse satellite data at irregular time step. We propose to combine various methods and algorithms to estimate the parameters efficiently. Indeed, simulation results indicate that the method of moment leads to a computationally efficient and numerically robust estimation procedure suitable for initializing the EM algorithm. A standard numerical optimization procedure is then used in the vicinity of the maximum of the likelihood function identified by the EM algorithm. It permits to accelerate the convergence of the EM algorithm with the extra benefit of giving as output an estimate of the information matrix which provide an estimate of the variance of the estimates.

This paper focus on SST data from the METOP satellite and the various results given in this paper indicate that the model is appropriate for describing some important properties of this data set such as the temporal structure and the measurement errors. Comparison with buoy data indicates that there is work to be done in order to estimate realistic SST conditions from METOP data. Nevertheless, we think that the state-space formulation adopted in this work is an appropriate method. In order to reconstruct realistic SST maps, we plan to extend the formulation in space and time to handle SST data from various satellites with their own accuracies and space-

time resolutions. Indeed, using such formulation has several benefits. First, it allows modelling flexibility. For example, non-linear dynamics, which incorporate the effects of advection and diffusion (see [20] and references therein) or non-linear evolution in the atmospheric variability can be considered. We also plan to investigate more elaborated measurement equations and include covariates to model the changing biases and variances of the different satellites (see e.g. [25]). Then, the Markovian structure of the model leads to efficient methods for the statistical inference. In particular, it allows to compute the maximum likelihood estimates and it is shown that these estimates are more efficient than the ones obtained using the method of moment commonly used in geostatitics with kriging. The Kalman recursions used to compute the smoothing probabilities take also benefit of the Markovian properties of the model and permit to save computational time compared to space-time kriging where high dimensional linear systems need to be solved.

**Acknowledgements** We would like to thank the TAO Project Office of the National Oceanic and Atmospheric Administration/ Pacific Marine Environmental laboratory (NOAA/PMEL), the National Climatic Data Center (NCDC) and the Godae High Resolution Sea Surface Temperature Pilot Project (GHRSST-PP) for respectively providing in situ PIRATA, real-time OIV2 SST analysis and satellite METOP SST measurements. We are grateful to the Meteo-France Lannion Team, A. Bentamy and O. Talagrand for their expertise and valuable comments on this work.
## References

1. Bertino L, Evensen G, Wackernagel H (2003) Sequential data assimilation techniques in oceanography. Int Stat Rev 71: 223-241

2. Cappé O, Moulines E, Rydén T (2005) Inference in hidden Markov models. Springer Verlag, New York

3. Cressie NAC (1993) Statistics for spatial data. Wiley, New York

4. Dempster AP, Laird NM, Rubin DB (1977) Maximum likelihood from incomplete data via the EM algorithm. J R Stat Soc B (Methodol) 39(1):1-38

5. Deng L, Shen X (1997) Maximum likelihood in statistical estimation of dynamic systems: decomposition algorithm and simulation results. Signal Process 57(1):65-79

6. Digalakis V, Rohlicek JR, Ostendorf M (1993) MI estimation of a stochastic linear system with the EM algorithm and its application to speech recognition 1(4):431-442

7. Donlon CJ, Minnett PJ, Gentemann C, Nightingale TJ, Barton IJ, Ward B, Murray MJ (2002) Toward improved validation of satellite sea surface skin temperature measurements for climate research. J Clim 15:353-369

8. Durbin J, Koopman SJ (2001) Time series analysis by state space methods. Oxford University Press, Oxford

9. Foltz GR, McPhaden MJ (2008) Impact of Saharan dust on tropical North Atlantic SST. J Clim 21:5048-5060

10. Frankignoul C, Hasselmann K (1977) Stochastic climate models. II Application to sea-surface temperature anomalies and thermocline variability. Tellus 29:289-305

11. Ghil M, Malanotte-Rizzoli P (1991) Data assimilation in meteorology and oceanography. Adv Geophys 33:141-266

12. Grodsky SA, Carton JA, Provost C, Servain J, Lorenzetti JA, McPhaden MJ (2005) Tropical instability waves at 0 N, 23 W in the Atlantic: a case study using Pilot Research Moored Array in the Tropical Atlantic (PIRATA) mooring data. J Geophys Res 110(C8):C08010

13. Heemink AW, Segers AJ (2002) Modeling and prediction of environmental data in space and time using Kalman filtering. Stoch Environ Res Risk Assess 16(3):225-240

14. Ide K, Courtier P, Ghil M, Lorenc AC (1997) Unified notation for data assimilation: operational, sequential and variational. Practice 75(1B):181-189

15. Jeff Wu CF (1983) On the convergence properties of the EM algorithm. Ann Stat 11(1):95-103

16. Jones RH, Boadi-Boateng F (1991) Unequally spaced longitudinal data with AR(1) serial correlation. Biometrics 47(1):161-175

17. Kaplan A, Cane MA, Kushnir Y, Clement AC, Blumenthal MB, Rajagopalan B (1998) Analyses of global sea surface temperature 1856-1991. J Geophys Res 103(18):567-589

18. Le Borgne P, Legendre G, Marsouin A (2007) Operational SST retrieval from MetOp/AVHRR. In: Proceedings of 2007 EU-METSAT conf., Amsterdam, The Netherlands

19. Penland C (1996) A stochastic model of Indo-Pacific sea surface temperature anomalies. Phys D Nonlinear Phenom 98(2-4):534-558

20. Piterbarg LI, Ostrovskii AG (1997) Advection and diffusion in random media: implications for sea surface temperature anomalies. Kluwer Academic Publishers, Dordrecht

21. Reynolds RW (1978) Sea surface temperature anomalies in the North Pacific Ocean. Tellus 30:97-103

22. Reynolds RW, Smith TM, Liu C, Chelton DB, Casey KS, Schlax MG (2007) Daily high-resolution-blended analyses for sea surface temperature. J Clim 20(22):5473-5496

23. Servain J, Busalacchi AJ, McPhaden MJ, Moura AD, Reverdin G, Vianna M, Zebiak SE (1998) A Pilot Research Moored Array in the Tropical Atlantic (PIRATA). Bull Am Meteorol Soc 79:2019-2032

24. Shumway RH, Stoffer DS (1982) An approach to time series smoothing and forecasting using the EM algorithm. J Time Ser Anal 3(4):253-264

25. Tandeo P, Autret E, Piolle JF, Tournadre J, Ailliot P (2009) A multivariate regression approach to adjust Aatsr sea surface temperature to in situ measurements. Geosci Remote Sens Lett IEEE 6(1):8-12

26. Wikle CK, Cressie N (1999) A dimension-reduced approach to space-time Kalman filtering. Biometrika 86:815-829