# NP – Nonlinear Processes in Geosciences

## NP1 – Mathematics for Planet Earth

### NP1.1 – Mathematics of Planet Earth

**Georg Gottwald**and Caroline Wormell

The long-term average response of observables of chaotic systems to dynamical perturbations can often be predicted using linear response theory, but not all chaotic systems possess a linear response. Macroscopic observables of complex dissipative chaotic systems, however, are widely assumed to have a linear response even if the microscopic variables do not, but the mechanism for this is not well-understood.

We present a comprehensive picture for the linear response of macroscopic observables in high-dimensional coupled deterministic dynamical systems, where the coupling is via a mean field and the microscopic subsystems may or may not obey linear response theory. We derive stochastic reductions of the dynamics of these observables from statistics of the microscopic system, and provide conditions for linear response theory to hold in finite dimensional systems and in the thermodynamic limit. In particular, we show that for large systems of finite size, linear response is induced via self-generated noise.

We present examples in the thermodynamic limit where the macroscopic observable satisfies LRT, although the microscopic subsystems individually violate LRT, as well a converse example where the macroscopic observable does not satisfy LRT despite all microscopic subsystems satisfying LRT when uncoupled. This latter, maybe surprising, example is associated with emergent non-trivial dynamics of the macroscopic observable. We provide numerical evidence for our results on linear response as well as some analytical intuition.

**Francesco Ragone**and Freddy Bouchet

Extreme events are a major topic of interest in climate science. Studying rare extreme events with complex numerical climate models is computationally challenging, since in principle very long simulations are needed to sample a sufficient number of events to provide a reliable statistics. This problem can be tackled using rare event algorithms, numerical tools developed in the past decades in mathematics and statistical physics, dedicated to the reduction of the computational effort required to sample rare events in dynamical systems. Typically they are designed as genetic algorithms, in which a set of suppression and cloning rules are applied to an ensemble simulation in order to focus the computational effort on the trajectories leading to the events of interest. Recently we showed the great potential of rare event algorithms for climate modelling, applying a rare event algorithm to study extremes of European surface temperature in Plasim, an intermediate complexity model, in absence of external time dependent forcings (no seasonal and daily cycles). Here we go beyond these limitations, studying extreme heat waves and warm summers in the Northern extratropics in fully realistic conditions including daily and seasonal cycles, both in Plasim and in the state of the art Earth system model CESM. We show how the algorithm allows to sample extreme events characterised by persistency on different time scales, discussing links with large deviation theory. We show how one can characterise the statistics of heat waves and warm summers with extremely large return times, with computational costs orders of magnitude smaller than with direct sampling, and reach ultra rare events that would have been impossible to observe otherwise. We analyse the emergence of teleconnection patterns during the extreme events and their relation to the dynamics of planetary waves. Finally we discuss how these results open the way to the systematic application of these techniques to a vast range of applicative and theoretical studies.

**Jakob Runge**

Detecting causal relationships from observational time series datasets is a key problem in better understanding the complex dynamical system Earth. Recent methodological advances have addressed major challenges such as high-dimensionality and nonlinearity, e.g., PCMCI (Runge et al. Sci. Adv. 2019), but many more remain. In this talk I will give an overview of challenges and methods and present a novel algorithm to identify causal directions among contemporaneous (or instantaneous) relationships. Such contemporaneous relations frequently appear when time series are aggregated (e.g., at a monthly resolution). Then approaches such as Granger Causality and PCMCI fail because they currently only address time-lagged causal relations.

We present extensive numerical examples and results on the causal relations among major climate modes of variability. The work overcomes a major drawback of current causal discovery methods and opens up entirely new possibilities to discover causal relations from time series in climate research and other fields in geosciences.

Runge et al., Detecting and quantifying causal associations in large nonlinear time series datasets, Science Advances eeaau4996 (2019).

**Lea Oljaca**, Jochen Broecker, and Tobias Kuna

Data assimilation is a term used to describe efforts to improve our knowledge

of a system by combining incomplete observations with imperfect models.

This is more generally known as filtering, which is ’optimal’ estimation of

the state of a system as it evolves over time, in the mean square sense. In

a Bayesian framework, the optimal filter is therefore naturally a sequence of

conditional probabilities of a signal given the observations and can be up-

dated recursively with new observations with Bayes’ formula. When, the

dynamics and observations errors are linear, this is equivalent to the Kalman

filter. In the nonlinear case, deriving an explicit form for the posterior dis-

tribution is in general not possible.

One of the important difficulties with applying the nonlinear filter in practice

is that the initial condition, the prior, is required to initialise the filtering.

However we are unlikely to know the correct initial distribution accurately

or at all. A filter is called stable if it is insensitive with respect to the

prior, that is, it converges to the same distribution, regardless of the initial

condition.

A body of work exists showing stability of the filter which rely on the stochas-

ticity of the underlying dynamics. In contrast, we show stability of the op-

timal filter for a class of nonlinear and deterministic dynamical systems and

our result relies on the intrinsic chaotic properties of the dynamics. We build

on the considerable knowledge that exists on the existence of SRB measures

in uniformly hyperbolic dynamical systems and we view the conditional prob-

abilities as SRB measures ‘conditional on the observation’ which are shown

to be absolutely continuous along the unstable manifold. This is in line with

the result of Bouquet, Carrassi et al [1] regarding data assimilation in the

“unstable subspace”, where they show stability of the filter if the unstable

and neutral subspaces are uniformly observed.

[1] M. Bocquet et al. “Degenerate Kalman Filter Error Covariances and

Their Convergence onto the Unstable Subspace”. In: SIAM/ASA Jour-

nal on Uncertainty Quantification 5.1 (2017), pp. 304–333. url: https:

//doi.org/10.1137/16M1068712.

**Julien Brajard**, Alberto Carrassi, Marc Bocquet, and Laurent Bertino

Can we build a machine learning parametrization in a numerical model using sparse and noisy observations?

In recent years, machine learning (ML) has been proposed to devise data-driven parametrizations of unresolved processes in dynamical numerical models. In most of the cases, ML is trained by coarse-graining high-resolution simulations to provide a dense, unnoisy target state (or even the tendency of the model).

Our goal is to go beyond the use of high-resolution simulations and train ML-based parametrization using direct data. Furthermore, we intentionally place ourselves in the realistic scenario of noisy and sparse observations.

The algorithm proposed in this work derives from the algorithm presented by the same authors in https://arxiv.org/abs/2001.01520.The principle is to first apply data assimilation (DA) techniques to estimate the full state of the system from a non-parametrized model, referred hereafter as the physical model. The parametrization term to be estimated is viewed as a model error in the DA system. In a second step, ML is used to define the parametrization, e.g., a predictor of the model error given the state of the system. Finally, the ML system is incorporated within the physical model to produce a hybrid model, combining a physical core with a ML-based parametrization.

The approach is applied to dynamical systems from low to intermediate complexity. The DA component of the proposed approach relies on an ensemble Kalman filter/smoother while the parametrization is represented by a convolutional neural network.

We show that the hybrid model yields better performance than the physical model in terms of both short-term (forecast skill) and long-term (power spectrum, Lyapunov exponents) properties. Sensitivity to the noise and density of observation is also assessed.

**Giulia Carigi**, Jochen Bröcker, and Tobias Kuna

In the Climate Sciences, there is great interest in understanding the long term average behaviour of the climate system. In the context of climate models, this behaviour can be expressed intrinsically by concepts from the theory of dynamical systems such as attractors and invariant measures. In particular to ensure long term statistics of the model are well defined from a mathematical perspective, the model needs to admit a unique ergodic invariant probability measure.

In this work we consider a classic two layer quasi geostrophic model, with the upper layer perturbed by additive noise, white in time and coloured in space, in order to account for random forcing, for instance through wind shear. Existence and uniqueness of an ergodic invariant measure is established using a technique from stochastic analysis called asymptotic coupling as described in [1]. The main difficulty in the proof is to show that two copies of the system that are driven by the same noise realisation can be synchronised through a coupling. This coupling has to be finite dimensional and act only on the upper layer.

Our results will be a key step to allow rigorous investigation of the response theory for the system under concern.

[1] Glatt-Holtz, N., Mattingly, J.C. & Richards, G. J Stat Phys (2017) 166: 618.

**Michael Ghil**, Gisela D. Charó, Denisse Sciamarella, and Mickael D. Chekroun

Chekroun et al. (*Physica D*, **240**, 2011) studied the global random attractor associated with the Lorenz (1963) model driven by multiplicative noise; they dubbed this time-evolving attractor LORA for short. The present talk examines the topological structure of the snapshots that approximate LORA’s evolution in time.

Sciamarella & Mindlin (*Phys. Rev. Lett., ***82**, 1999;* Phys. Rev. E*, **64**, 2001) introduced the methodology of Branched Manifold Analysis through Homologies (BraMAH) to the study of chaotic flows. Here, the BraMAH methodology is extended for the first time, to the best of our knowledge, from deterministically chaotic flows to nonlinear noise-driven systems.

The BraMAH algorithm starts from a cloud of points given by a large number of orbits and it builds a rough skeleton of the underlying branched manifold on which the points lie. This construction is achieved by local approximations of the manifold that use Euclidean closed sets; the nature of these sets depends on their topological dimension, e.g., segments or disks. The skeleton is mathematically expressed as a complex of cells, whose algebraic topology is analyzed by computing its homology groups.

The analysis is performed for a fixed realization of the driving noise at different time instants. We show that the topology of LORA changes in time and that it is quite distinct from the time-independent one of the classical Lorenz (1963) convection model, for the same values of the parameters. Topological tipping points are also studied by varying the parameter values from the classical ones.

**Manita Chouksey**

Geophysical flows such as the atmosphere and the ocean are characterized by rotation and stratification, which together give rise to two dominant motions: the slow balanced and the fast unbalanced motions. The interaction between the balanced and unbalanced motions and the energy transfers between them impact the energy and momentum cycle of the flow, and is therefore crucial to understand the underlying energetics of the atmosphere and the ocean. Balanced motions, for instance mesoscale eddies, can transfer their energy to unbalanced motions, such as internal gravity waves, by spontaneous loss of balance amongst other processes. The exact mechanism of wave generation, however, remain less understood and is hindered to an extent by the challenge of separating the flow field into balanced and unbalanced motions.

This separation is achieved using two different balancing procedures in an identical model setup and assess the differences in the obtained balanced state and the resultant energy transfer to unbalanced motions. The first procedure we implement is a non-linear initialisation procedure based on Machenhauer (1977) but extended to higher orders in Rossby number. The second procedure implemented is the optimal potential vorticity balance to achieve the balanced state. The results show that the numerics of the model affect the obtained balanced state from the two procedures, and thus the residual signal which we interpret as the unbalanced motions, i.e. internal gravity waves. A further complication is the presence of slaved modes, which appear along the unbalanced motions but are tied to the balanced motions, for which we need to extend the separation to higher orders in Rossby number. Further, we assess the energy transfers between balanced and unbalanced motions in experiments with different Rossby numbers and for different orders in Rossby number. We find that it is crucial to consider the effect of the numerics in models and make a suitable choice of the balancing procedure, as well as diagnose the unbalanced motions at higher orders to precisely detect the unbalanced wave signal.

**Maximilian Gelbrecht**, Jürgen Kurths, and Frank Hellmann

Many high-dimensional complex systems such as climate models exhibit an enormously complex landscape of possible asymptotic state. On most occasions these are challenging to analyse with traditional bifurcation analysis methods. Often, one is also more broadly interested in classes of asymptotic states. Here, we present a novel numerical approach prepared for analysing such high-dimensional multistable complex systems: Monte Carlo Basin Bifurcation Analysis (MCBB). Based on random sampling and clustering methods, we identify the type of dynamic regimes with the largest basins of attraction and track how the volume of these basins change with the system parameters. In order to due this suitable, easy to compute, statistics of trajectories with randomly generated initial conditions and parameters are clustered by an algorithm such as DBSCAN. Due to the modular and flexible nature of the method, it has a wide range of possible applications. While initially oscillator networks were one of the main applications of this methods, here we present an analysis of a simple conceptual climate model setup up by coupling an energy balance model to the Lorenz96 system. The method is available to use as a package for the Julia language.

**Shaun Lovejoy**, Lenin Del Rio Amador, and Roman Procyk

The conventional 1-D energy balance equation (EBE) has no vertical coordinate so that radiative imbalances between the earth and outer space are redirected in the horizontal in an ad hoc manner. We retain the basic EBE but add a vertical coordinate so that the imbalances drive the system by imposing heat fluxes through the surface. While this is theoretically correct, it leads to (apparently) difficult mixed boundary conditions. However, using Babenko’s method, we directly obtain simple analytic equations for (2D) surface temperature anomalies T_{s}(x,t): the Half-order Energy Balance Equation (HEBE) and the Generalized HEBE, (GHEBE) [Lovejoy, 2019a]. The HEBE anomaly equation only depends on the local climate sensitivities and relaxation times. We analytically solve the HEBE and GHEBE for T_{s}(x,t) and provide evidence that the HEBE applies at scales >≈10km. We also calculate very long time diffusive transport dominated climate states as well as space-time statistics including the cross-correlation matrix needed for empirical orthogonal functions.

The HEBE is the special H = 1/2 case of the Fractional EBE (FEBE) [Lovejoy, 2019b], [Lovejoy, 2019c] and has a long (power law) memory up to its relaxation time t. Several papers have empirically estimated H ≈ 0.5, as well as t ≈ 4 years, whereas the classical zero-dimensional EBE has H = 1 and t ≈ 4 years. The former values permit accurate macroweather forecasts and low uncertainty climate projections; this suggests that the HEBE could apply to time scales as short as a month. Future generalizations include albedo-temperature feedbacks and more realistic treatments of past and future climate states.

**References**

Lovejoy, S., The half-order energy balance equation, J. Geophys. Res. (Atmos.), (submitted, Nov. 2019), 2019a.

Lovejoy, S., Weather, Macroweather and Climate: our random yet predictable atmosphere, 334 pp., Oxford U. Press, 2019b.

Lovejoy, S., Fractional Relaxation noises, motions and the stochastic fractional relxation equation Nonlinear Proc. in Geophys. Disc., https://doi.org/10.5194/npg-2019-39, 2019c.

**Michael Coughlan**, Ian Hewitt, Sam Howison, and Andrew Wells

Arctic sea ice forms a thin but significant layer at the ocean surface, mediating key climate feedbacks. During summer, surface melting produces considerable volumes of water, which collect on the ice surface in ponds. These ponds have long been suggested as a contributing factor to the discrepancy between observed and predicted sea ice extent. When viewed at large scales ponds have a complicated, approximately fractal geometry and vary in area from tens to thousands of square meters. Increases in pond depth and area lead to further increases in heat absorption and overall melting, contributing to the ice-albedo feedback.

Previous modelling work has focussed either on the physics of individual ponds or on the statistical behaviour of systems of ponds. We present a physically-based network model for systems of ponds which accounts for both the individual and collective behaviour of ponds. Each pond initially occupies a distinct catchment basin and evolves according to a mass-conserving differential equation representing the melting dynamics for bare and water-covered ice. Ponds can later connect together to form a network with fluxes of water between catchment areas, constrained by the ice topography and pond water levels.

We use the model to explore how the evolution of pond area and hence melting depends on the governing parameters, and to explore how the connections between ponds develop over the melt season. Comparisons with observations are made to demonstrate the ways in which the model qualitatively replicates properties of pond systems, including fractal dimension of pond areas and two distinct regimes of pond complexity that are observed during their development cycle. Different perimeter-area relationships exist for ponds in the two regimes. The model replicates these relationships and exhibits a percolation transition around the transition between these regimes, a facet of pond behaviour suggested by previous studies. Our results reinforce the findings of these studies on percolation thresholds in pond systems and further allow us to constrain pond coverage at this threshold - an important quantity in measuring the scale and effects of the ice-albedo feedback.

**Michael Berhanu**, Raphael Dubourg, Arthur Walbecq, Cyril Ozouf, Adrien Guerin, Julien Derr, and Sylvain Courrech du Pont

Erosion by dissolution is a decisive process shaping small-scale landscape morphology [1]. For fast dissolving minerals, the erosion rate is controlled by the solute transport [2] and characteristic erosion patterns can appear due to hydrodynamics mechanisms. Among the diversity of the dissolution patterns, the scallops are small depressions in a dissolving wall, appearing as cups with sharp edges. Their size varies from few millimeters to around ten centimeters. The scallops occur typically as the final steady form of ripple patterns created by the action of a turbulent flow on a dissolving surface [3,4]. Moreover, very similar shapes are also met, without imposed external flow, when the fluid motion results from the solutal convection induced by the dissolution [2,5,6]. Finally, scallop patterns resulting from similar mechanisms appear also on ice surfaces by melting in presence of a turbulent flow [7] or a convection flow [6].

Using three-dimensional surface reconstruction, we characterize quantitatively the scallop patterns mainly for experimental samples patterned by solutal convection. The temporal evolution of the scallop shape, of their spatial distribution and of the induced roughness are specifically investigated, in order to determine mechanisms explaining the generic aspects of dissolution patterns.

[1] P. Meakin and B. Jamtveit, Geological pattern formation by growth and dissolution in aqueous systems, **Proc. R. Soc. A 466** 659-694 (2010)

[2] J. Philippi, M. Berhanu, J. Derr and S. Courrech du Pont, Solutal convection induced by dissolution, **Phys. Rev. Fluids, 4,** 103801 (2019)

[3] P.N. Blumberg and R.L. Curl, Experimental and theoretical studies of dissolution roughness, **J. Fluid Mech. 65**, 735 (1974)

[4] P. Claudin, O. Durán and B. Andreotti, Dissolution instability and roughening transition, **J. Fluid Mech. 832**, R2 (1974)

[5] T.S. Sullivan, Y. Liu and R. E. Ecke, Turbulent solutal convection and surface patterning in solid dissolution, **Phys. Rev. E 54**, (1) 486, (1996)

[6] C. Cohen, M. Berhanu, J. Derr and S. Courrech du Pont, Erosion patterns on dissolving and melting bodies (2015 Gallery of Fluid motion), **Phys. Rev. Fluids, 1,** 050508 (2016)

[7] M. Bushuk, D. M. Holland, T. P. Stanton, A. Stern and C. Gray. Ice scallops: a laboratory investigation of the Ice-water interface, **J. Fluid Mech. 873**, 942 (2019)

**Pauline Delorme**, Giles Wiggs, Matthew Baddock, Joanna Nield, James Best, Kenneth Christensen, Nathaniel Bristow, Andrew Valdez, and Philippe Claudin

Early-stage aeolian bedforms develop into sand dunes through complex interactions between flow, sediment transport and surface topography. Depending on the specific environmental and wind conditions the mechanisms of dune formation, and ultimately the shape of the nascent dunes, may differ. Here, we investigate the formation of a proto-dune-field, located in the Great Sand Dunes National Park ( Colorado, USA), using a three dimensional linear stability analysis.

We use in-situ measurements of wind and sediment transport, collected during a one-month field campaign, as part of a linear stability analysis to predict the orientation and wavelength of the proto-dunes.

We find that the output of the linear stability analysis compares well to high-resolution Digital Elevation Models measured using terrestrial laser scanning. Our findings suggest that the bed instability mechanism is a good predictor of proto-dune development on sandy surfaces with a bimodal wind regime.

**Milad Hooshyar**, Sara Bonetti, Arvind Singh, Efi Foufoula-Georgiou, and Amilcare Porporato

The channelization cascade observed in terrestrial landscapes describes the progressive formation of large channels from smaller ones starting from diffusion-dominated hillslopes. This behavior is reminiscent of other non-equilibrium complex systems, particularly fluids turbulence, where larger vortices break down into smaller ones until viscous dissipation dominates. Based on this analogy, we show that topographic surfaces emerging between parallel zero-elevation boundaries present a logarithmic scaling in the mean-elevation profile, which resembles the well-known logarithmic velocity profile in wall-bounded turbulence. Within this region of elevation fluctuation, the power spectrum exhibits a power-law decay resembling the Kolmogorov -5/3 scaling of turbulence. We also demonstrate that similar scaling behaviors emerge in surfaces from a laboratory experiment, natural basins, and constructed following optimality principles. In general, we show that the steady-state solutions of the governing equations of landscape evolution are the stationary surfaces of a functional defined as the average domain elevation. Depending on the exponent of the specific drainage area in the erosion term (m), the steady-state surfaces are local minimum (m<1) or maximum (m>1) of the average domain elevation.

**Manuel Santos Gutiérrez**and Valerio Lucarini

Dynamical systems are often subject to forcing or changes in their governing parameters and it is of interest to study

how this affects their statistical properties. A prominent real-life example of this class of problems is the investigation

of climate response to perturbations. In this respect, it is crucial to determine what the linear response of a system is

as a quantification of sensitivity. Alongside previous work, here we use the transfer operator formalism to study the

response and sensitivity of a dynamical system undergoing perturbations. By projecting the transfer operator onto a

suitable finite dimensional vector space, one is able to obtain matrix representations which determine finite Markov

processes. Further, using perturbation theory for Markov matrices, it is possible to determine the linear and nonlinear

response of the system given a prescribed forcing. Here, we suggest a methodology which puts the scope on the

evolution law of densities (the Liouville/Fokker-Planck equation), allowing to effectively calculate the sensitivity and

response of two representative dynamical systems.

**Valerio Lembo**, Valerio Lucarini, and Francesco Ragone

Global Climate Models are key tools for predicting the future response of the climate system to a variety of natural and anthropogenic forcings. Typically, an ensemble of simulations is performed considering a scenario of forcing, in order to analyse the response of the climate system to the specific forcing signal. Given that the the climate response spans a very large range of timescales, such a strategy often requires a dramatic amount of computational resources. In this paper we show how to use statistical mechanics to construct operators able to flexibly predict climate change for a variety of climatic variables of interest, going beyond the limitation of having to consider specific time patterns of forcing. We perform our study on a fully coupled GCM - MPI-ESM v.1.2 - and for the first time we prove the effectiveness of response theory in predicting future climate response to CO_{2} increase on a vast range of temporal scales. We specifically treat atmospheric (surface temperature) and oceanic variables (strength of the Atlantic Meridional Overturning Circulation and of the Antarctic Circumpolar Current), as well as the global ocean heat uptake.

**Hege-Beate Fredriksen**

We investigate a new algorithm for estimating time-evolving global forcing in climate models. The method is an extension of a previous method by Forster et al. (2013), but here we also allow for a globally nonlinear feedback. We assume the nonlinearity of this global feedback can be explained as a time-scale dependence, associated with linear temperature responses to the forcing on different time scales, as in Proistosescu and Huybers (2017). With this method we obtain stronger forcing estimates than previously believed for the representative concentration pathway experiments in CMIP5 models. The reason for the higher future forcing is that the global feedback has a higher magnitude at the smaller time scales than at the longer time scales – this is closely related to provided arguments for the equilibrium climate sensitivity showing changes with time.

We find also that the linear temperature response to our new forcing predicts the modelled response quite well, although the response is a little overestimated for some CMIP5 models. Finally, we discuss the assumptions made in our study, and consistency of our assumptions with other leading hypotheses for why the global feedback is nonlinear.

References:

Forster, P. M., T. Andrews, P. Good, J. M. Gregory, L. S. Jackson, and M. Zelinka (2013), Evaluating adjusted forcing and model spread for historical and future scenarios in the cmip5 generation of climate models, Journal of Geophysical Research, 118, 1139–1150, doi:10.1002/jgrd.50174.

Proistosescu, C., and P. J. Huybers (2017), Slow climate mode reconciles historical and model-based estimates of climate sensitivity, Sci. Adv., 3, e1602, 821, doi:10.1126/sciadv.1602821

**Chiara Cecilia Maiocchi**, Valerio Lucarini, Andrey Gritsun, and Grigorios Pavliotis

Unstable periodic orbits (UPOs) have been proved to be a relevant mathematical tool in the study of Climate Science. In a recent paper Lucarini and Gritsun [1] provided an alternative approach for understanding the properties of the atmosphere. Climate can be interpreted as a non-equilibrium steady state system and, as such, statistical mechanics can provide us with tools for its study.

UPOs decomposition plays a relevant role in the study of chaotic dynamical systems. In fact, UPOs densely populate the attractor of a chaotic system, and can therefore be thought as building blocks to construct the dynamic of the system itself. Since they are dense in the attractor, it is always possible to find a UPO arbitrarily near to a chaotic trajectory: the trajectory will remain close to the UPO, but it will never follow it indefinitely, because of its instability. Loosely speaking, a chaotic trajectory is repelled between neighbourhoods of different UPOs and can thus be approximated in terms of these periodic orbits. The characteristics of the system can then be reconstructed from the full set of periodic orbits in this fashion.

The sampling of UPOs is therefore a relevant problem for describing chaotic dynamical systems and can represent an interesting topic for the study of Climate Science. In this work we address this problem and present an algorithm to numerically extract UPOs from the attractor of a simple Climate Model such as Lorenz-63.

[1] V. Lucarini and A. Gritsun, “A new mathematical framework for atmospheric blocking events,” Climate Dynamics, vol. 54, pp. 575–598, Jan 2020.

**Michel Crucifix**, Dmitri Alexandrov, irina Bashkirtseva, and Lev Ryashko

Glacial-interglacial cycles are global climatic changes which have characterised the last 3 million years. The eight latest

glacial-interglacial cycles represent changes in sea level over 100 m, and their average duration was around 100 000 years. There is a

long tradition of modelling glacial-interglacial cycles with low-order dynamical systems. In one view, the cyclic phenomenon is caused by

non-linear interactions between components of the climate system: The dynamical system model which represents Earth dynamics has a limit cycle. In an another view, the variations in ice volume and ice sheet extent are caused by changes in Earth's orbit, possibly amplified by feedbacks.

This response and internal feedbacks need to be non-linear to explain the asymmetric character of glacial-interglacial cycles and their duration. A third view sees glacial-interglacial cycles as a limit cycle synchronised on the orbital forcing.

The purpose of the present contribution is to pay specific attention to the effects of stochastic forcing. Indeed, the trajectories

obtained in presence of noise are not necessarily noised-up versions of the deterministic trajectories. They may follow pathways which

have no analogue in the deterministic version of the model. Our purpose is to

demonstrate the mechanisms by which stochastic excitation may generate such large-scale oscillations and induce intermittency. To this end, we

consider a series of models previously introduced in the literature, starting by autonomous models with two variables, and then three

variables. The properties of stochastic trajectories are understood by reference to the bifurcation diagram, the vector field, and a

method called stochastic sensitivity analysis. We then introduce models accounting for the orbital forcing, and distinguish forced and

synchronised ice-age scenarios, and show again how noise may generate trajectories which have no immediate analogue in the determinstic model.

**Nikolas Porz**

The representation of cloud processes in weather and climate models is crucial for their feedback on atmospheric flows. Since there is no general macroscopic theory of clouds, the parameterization of clouds in corresponding simulation software depends fundamentally on the underlying modeling assumptions. We present a new model of intermediate complexity (a one-and-a-half moment scheme) for warm clouds, which is derived from physical principles. Our model consists of a system of differential-algebraic equations which allows for supersaturation and thus avoids the commonly used but somewhat outdated concept of so called 'saturation adjustment'. This is made possible by a non-Lipschitz right-hand side, which allows for nontrivial solutions. In a recent effort we have proved under mild assumptions on the external forcing that this system of equations has a unique physically consistent solution, i.e., a solution with a nonzero droplet population in the supersaturated regime. For the numerical solution of this system we have developed a semi-implicit integration scheme, with efficient solvers for the implicit parts. The model conserves air and water (if one accounts for the precipitation), and it comes with eight parameters that cannot be measured since they describe simplified processes, so they need to be fitted to the data. For further studies we implemented our cloud micro physics model into ICON, the weather forecast model operated by the German forecast center DWD.

**Robert Malte Polzin**, Annette Müller, Peter Nevir, Henning Rust, and Peter Koltai

The presented work contains an investigation of the stochastic aggregation of convective structures on different scales in the atmosphere. A

computational framework is applied that provides highly scalable identification of reduced Bayesian models. The deterministic large scale

flow variables are reduced into latent states, whereas the stochastic small scale convective structures are affiliated to these. The analysis of

the latent states in number and maximization reduction improves the understanding for the large scale forcing of convective processes. The

convective structures are determined by vertical velocities. Different variables of the large-scale flow, such as the convective available

potential energy, available moisture, vertical windshear and the Dynamic State Index (DSI), a diabaticity indicator, are investigated. Our approach

does not require a distributional assumption but works instead with a discretised and categorised state vector.

**Vera Melinda Galfi**, and Valerio Lucarini

We use large deviation theory to study persistent extreme events of temperature, like heat waves or cold spells. We consider the mid-latitudes of a simplified yet Earth-like general circulation model of the atmosphere and numerically estimate large deviation rate functions of near-surface temperature averages over different spatial scales. We find that, in order to represent persistent extreme events based on large deviation theory, one has to look at temporal averages of spatially averaged observables. The spatial averaging scale is crucial, and has to correspond with the scale of the event of interest. Accordingly, the computed rate functions indicate substantially different statistical properties of temperature averages over intermediate spatial scales (larger, but still of the order of the typical scale), as compared to the ones related to any other scale. Thus, heat waves (or cold spells) can be interpreted as large deviations of temperature averaged over intermediate spatial scales. Furthermore, we find universal characteristics of rate functions, based on the equivalence of temporal, spatial, and spatio-temporal rate functions if we perform a re-normalisation by the integrated auto-correlation.

**Tobias Kuna**, Valerio Lucarini, Davide Faranda, Jerouen Wouters, and Viviane Baladi

Extremes are related to high impact and serious hazard events and hence their study and prediction have been and continue to be highly relevant for all kind of applications in geoscience and beyond. Extreme value theory is promising to be able to predict them reliably and robustly. In the last fifteen years the classical extreme value theory for stochastic processes has been extended to dynamical systems and has been related to properties of physical measure (statistical properties of the system), return and hitting times. We will review what one can say for highly dimensional perfectly chaotic systems. We will concentrate on relations between the index of the extreme distribution and invariants of the underlying dynamical system which are stable, in the sense that they will continuously depend on changing parameters in the dynamics. Furthermore, we explore whether there exists a response theory for extremes, that is, whether the change of extremes can be quantitatilvely expressed in terms of changing parameters.

**Paul Bowen**

Approximations in the moist thermodynamics of atmospheric/weather models are often inconsistent. Different parts of numerical models may handle the thermodynamics in different ways, or the approximations may disagree with the laws of thermodynamics. In order to address these problems, we may derive all relevant thermodynamic quantities from a defined thermodynamic potential; approximations are then instead made to the potential itself — this guarantees self-consistency. This concept is viable for vapor and liquid water mixtures in a moist atmospheric system using the Gibbs function but on extension to include the ice phase an ambiguity presents itself at the triple-point. In order to resolve this the energy function must be utilised instead; constrained maximisation methods can then be used on the entropy in order to solve the system equilibrium state. Once this is done however, a further extension is necessary for atmospheric systems. In the Earth’s atmosphere many important non-equilibrium processes take place; for example, freezing of super-cooled water, evaporation, and precipitation. To fully capture these processes the equilibrium method must be reformulated to involve finite rates of approach towards equilibrium. This may be done using various principles of non-equilibrium thermodynamics, principally Onsager reciprocal relations. A numerical scheme may then be implemented which models competing finite processes in a moist thermodynamic system.

**Charline Ragon**, Valerio Lembo, Valerio Lucarini, Jérôme Kasparian, and Maura Brunetti

The climate system can be seen as a thermal engine that generates entropy by irreversible processes and achieves a steady state by redistributing the input solar energy among its different components (ocean, atmosphere, etc) and by balancing the energy, water mass and entropy budgets over all the spatial scales. Biases in modern climate models are generally related to the fact that their statistical properties are not well represented, giving rise to imperfect closures of the energy cycle. Thus, a proper measurement of the efficiency of the thermal engine in each climate model is needed. Moreover, possible steady states (attractors) that can be approached at climate tipping-points are characterised by different feedbacks becoming dominant in the thermal engine.

We apply the Thermodynamic Diagnostic Tool (*TheDiaTo*) [1] to the attractors recently obtained using the MIT general circulation model (*MITgcm*) in a coupled aquaplanet [2], a planet where the ocean covers the entire globe. Such coupled aquaplanets, where nonlinear interactions between atmosphere, ocean and sea ice are fully taken into account, provide a relevant framework to understand the role of the main feedbacks at play in the climate system. Five attractors have been found, ranging from snowball (where ice covers the entire planet) to hot state conditions (where ice completely disappears) [2].

Using *TheDiaTo*, we analyse the five climate attractors by estimating: a) the energy budgets and meridional energy transports; b) the water mass and latent energy budgets and respective meridional transports; c) the Lorenz Energy Cycle; d) the material entropy production. We consider different coupled atmosphere-ocean-sea ice configurations and cloud parameterizations of *MITgcm* where the energy balance at the top of the atmosphere is progressively better closed in order to understand the occurrence of possible biases in the statistical properties of each attractor.

Our contribution will help clarify the thermodynamic differences in climate attractors and their stability to external perturbations that could shift the climate from a steady state to the other.

[1] Lembo V., Lunkeit F., Lucarini V., TheDiaTo (v1.0) – a new diagnostic tool for water, energy amd entropy budgets in climate models, Geosci. Model Dev. 12, 3805-3834 (2019)

[2] Brunetti M., Kasparian J., Vérard C., Co-existing climate attractors in a coupled aquaplanet, Climate Dynamics 53, 6293-6308 (2019)

**Stan Schymanski**and Martijn Westhoff

Thermodynamic optimality principles, such as maximum entropy production or maximum power extraction, hold a great promise to help explain self-organisation of various compartments of planet Earth, including the climate system, catchments and ecosystems. There is a growing number of examples for more or less successful use of these principles in earth system science, but a common systematic approach to the formulation of the relevant system boundaries, state variables and exchange fluxes has not yet emerged. Here we present a blueprint for the thermodynamically consistent formulation of box models and rigorous testing of optimality principles, in particular the maximum entropy production (MEP) and the maximum power (MP) principle. We investigate under what conditions these principles can be used to predict energy transfer coefficients across internal system boundaries and demonstrate that, contrary to common perception, these principles do not lead to similar predictions if energy and entropy balances are explicitly considered for the whole system and the defined sub-systems. We further highlight various pitfalls that may result in thermodynamically inconsistent models and potentially wrong conclusions about the implications of thermodynamic optimality principles.

The analysis is performed in an open source mathematical framework, using the notebook interface Jupyter, the programming language Python, Sympy and a newly developed package for Python, "Environmental Science using Symbolic Math" (ESSM, https://github.com/environmentalscience/essm). This ensures easy verifiability of the results and enables users to re-use and modify variable definitions, equations and mathematical solutions to suit their own thermodynamic problems.

**Bianca Kretz**, Willi Freeden, and Volker Michel

The aspect of poroelasticity is anywhere interesting where a solid material and a fluid come into play and have an effect on each other. This is the case in many applications and we want to focus on geothermics. It is useful to consider this aspect since the replacement of the water in the reservoir below the Earth's surface has an effect on the sorrounding material and vice versa. The underlying physical processes can be described by partial differential equations, called the quasistatic equations of poroelasticity (QEP). From a mathematical point of view, we have a set of three (for two space and one time dimension) partial differential equations with the unknowns u (displacement) and p (pore pressure) depending on the space and the time.

Our aim is to do a decomposition of the data given for u and p in order that we can see underlying structures in the different decomposition scales that cannot be seen in the whole data.

For this process, we need the fundamental solution tensor of the QEP (cf. [1],[5]).

That means we assume that we have given data for u and p (they can be obtained for example by a method of fundamental solutions, cf. [1]) and want to investigate a post-processing method to these data. Here we follow the basic approaches for the Laplace-, Helmholtz- and d'Alembert equation (cf. [2],[4],[6]) and the Cauchy-Navier equation as a tensor-valued ansatz (cf. [3]). That means we want to modify our elements of the fundamental solution tensor in such a way that we smooth the singularity concerning a parameter set τ=(τ_{x},τ_{t}).

With the help of these modified functions, we construct scaling functions which have to fulfil the properties of an approximate identity.

They are convolved with the given data to extract more details of u and p.

**References**

[1] M. Augustin: A method of fundamental solutions in poroelasticity to model the stress field in geothermal reservoirs, PhD Thesis, University of Kaiserslautern, 2015, Birkhäuser, New York, 2015.

[2] C. Blick, Multiscale potential methods in geothermal research: decorrelation reflected post-processing and locally based inversion, PhD Thesis, Geomathematics Group, Department of Mathematics, University of Kaiserslautern, 2015.

[3] C. Blick, S. Eberle, Multiscale density decorrelation by Cauchy-Navier wavelets, Int. J. Geomath. 10, 2019, article 24.

[4] C. Blick, W. Freeden, H. Nutz: Feature extraction of geological signatures by multiscale gravimetry. Int. J. Geomath. 8: 57-83, 2017.

[5] A.H.D. Cheng and E. Detournay: On singular integral equations and fundamental solutions of poroelasticity. Int. J. Solid. Struct. 35, 4521-4555, 1998.

[6] W. Freeden, C. Blick: Signal decorrelation by means of multiscale methods, World of Mining, 65(5):304-317, 2013.

**Matteo Bernard Bertagni**and Carlo Camporeale

The interactions between water and rocks create an extensive variety of marvelous patterns, which span on several classes of time and space scales. In this work, we provide a mathematical model for the formation of longitudinal erosive patterns commonly found in karst and alpine environments. The model couples the hydrodynamics of a laminar flow of water (Orr-Somerfield equation) to the concentration field of the eroded-rock chemistry. Results show that an instability of the plane rock wetted by the water film leads to a longitudinal channelization responsible for the pattern formation. The spatial scales predicted by the model span over different orders of magnitude depending on the flow intensity and this may explain why similar patterns of different sizes are observed in nature (millimetric microrills, centimetric rillenkarren, decametric solution runnels).

**Imogen Dell**

There exists a coupling mechanism between the troposphere and the stratosphere, which plays a fundamental role in weather and climate. The coupling is highly complex and rests upon radiative and chemical feedbacks, as well as dynamical coupling by Rossby waves. The troposphere acts as a source of Rossby waves which propagate upwards in to the stratosphere, affecting the zonal mean flow. Rossby waves are also likely to play a significant role in downward communication of information via reflection from the stratosphere in to the troposphere. A mechanism for this reflection could be from a so-called critical layer. A shear flow exhibits a critical layer where the phase speed equals the flow velocity, where viscous and nonlinear effects become important. A wave incident upon a critical layer may be absorbed, reflected or overreflected, whereby the amplitude of the reflected wave is larger than that of the incident wave. In the case of troposphere-stratosphere coupling, the concept of critical layer overreflection is key to understanding atmospheric instability.

Motivated by this, a mathematical framework for understanding the coupling will be presented together with an investigation in to the role of nonlinearity versus viscosity inside the critical layer.

**Gözde Özden**and Marcel Oliver

Consider the motion of a rotating fluid governed by the Boussinesq equations with full Coriolis parameter. This is contrary to the so-called ''traditional approximation'' in which the horizontal part of the Coriolis parameter is zero. The model is obtained using variational principle which depends on Lagrangian dynamics. The full Coriolis force is used since the horizontal component of the angular velocity has a crucial role in that it introduces a dependence on the direction of the geostrophic flow in the horizontal geostrophical plane. We aim that singularity near the equatorial region can be solved with this assumption. This gives a consistent balance relation for any latitude on the Earth. We follow the similar strategy to that Oliver and Vasylkevych (2016) for the system to derive the Euler-Poincaré equations. Firstly, the system is transformed into desired scale giving the differences with the other scales. We derive the balance model Lagrangian as called L1 model, R. Salmon, using Hamiltonian principles. Near identity transformation is applied to simplify the Hamiltonian. Whole calculations are done considering the smallness assumption of the Rossby number. Long term, we aim that results help to understand the global energy cycle with the goal of validity and improving climate models.

**Stephen Griffiths**

The Kelvin wave is perhaps the most important of the equatorially trapped waves in the terrestrial atmosphere and ocean, and plays a role in various phenomena such as tropical convection and El Nino. Theoretically, it can be understood from the linear dynamics of a stratified fluid on an equatorial β-plane, which, with simple assumptions about the disturbance structure, leads to wavelike solutions propagating along the equator, with exponential decay in latitude. However, when the simplest possible background flow is added (with uniform latitudinal shear), the Kelvin wave (but not the other equatorial waves) becomes unstable. This happens in an extremely unusual way: there is instability for arbitrarily small nondimensional shear *λ*, and the growth rate is proportional to exp(-1/λ^2) as λ → 0. This in contrast to most hydrodynamic instabilities, in which the growth rate typically scales as a positive power of λ-λ_{c} as the control parameter λ passes through a critical value λ_{c}.

This Kelvin wave instability has been established numerically by Natarov and Boyd, who also speculated as to the underlying mathematical cause by analysing a quantum harmonic oscillator perturbed by a potential with a remote pole. Here we show how the growth rate and full spatial structure of the Kelvin wave instability may be derived using matched asymptotic expansions applied to the (linear) equations of motion. This involves an adventure with confluent hypergeometric functions in the exponentially-decaying tails of the Kelvin waves, and a trick to reveal the exponentially small growth rate from a formulation that only uses regular perturbation expansions. Numerical verification of the analysis is also interesting and challenging, since special high-precision solutions of the governing ordinary differential equations are required even when the nondimensional shear is not that small (circa 0.5).

**Richard Blender**and Joscha Fregin

We consider recharge-discharge processes in a forced wave-mean flow interaction model and in a forced Rossby wave triad. Such processes are common in atmospheric dynamics and are typically modelled by nonlinear oscillators, for example for mid-latitude storms by Ambaum and Novak (2013) and for convective cycles by Yano and Plant (2012). A similar behaviour can be seen in the simulation of a forced wave number triad by Lynch (2009). Here we construct noncanonical Hamiltonian and Nambu representations in three-dimensional phase space for available and prescribed conservation laws during the recharge and discharge regimes. Divergence in phase space is modelled by a pre-factor. The approach allows the design of conservative and forced dynamical systems.

**Stuart Grieve**, Simon Mudd, Fiona Clubb, Michael Singer, Katerina Michaelides, and Shiuan-An Chen

The topology of fluvial networks has long been studied, with Horton's laws describing relationships between stream order, stream density, and stream length often cited as fundamental governing principles of drainage basin development. Building upon these principles, small scale studies have identified patterns of self-similarity in drainage networks in the continental USA, suggesting that to some extent, river networks self-organise in a scale invariant manner. More stringent measures of self-similarity have also been developed, which quantify the fractal nature of side branching structures in fluvial networks. Using such metrics, studies have identified similarities between leaf vein structures and fluvial networks, and have identified a potential climatic signature in North American river topology.

The appeal of such techniques over traditional methods of channel analysis using topographic data is that in self-similar networks, the precise location of channel heads is unimportant, allowing analysis to be performed at unprecedented scales, and in locations where data quality is limited. Here, we attempt to reconcile these two suites of techniques to understand the potential and limitations of network topology as an indicator of broader landscape dynamics. We achieve this through the analysis of fluvial networks extracted at a global scale from the Shuttle Radar Topography Mission dataset alongside other global earth observation data.

**Nathaniel Bristow**, James Best, Kenneth Christensen, Matthew Baddock, Giles Wiggs, Pauline Delorme, and Joanna Nield

Understanding the initiation of aeolian dunes poses significant challenges due to the strong couplings between turbulent fluid flow, sediment transport, and bedform morphology. While much is known concerning the dynamics of more mature bedforms, open questions remain as to how protodunes are formed, as well as the mechanisms by which they continue to evolve. The structure of the turbulent flow field drives the mobilization or deposition of sediment, thus controlling the initial formation of sand patches, yet is also strongly influenced itself by local conditions such as surface roughness and moisture. Furthermore, an additional feedback on the flow and transport is exerted by the sand patches themselves once they begin to form.

As protodunes begin to develop from this initial deposition, their morphologies possess unique characteristics involving a reverse asymmetry of the stoss and lee sides, wherein the crest begins upstream, close to the toe, and gradually shifts downstream toward the "regular" asymmetric profile exhibited by more mature dunes. However, these early stages of development also involve very gentle slopes and low profiles which make field measurements of the associated flow particularly challenging.

The current research effort involves a combination of field measurements, documenting the initiation and morphological development of sand patches and protodunes, in concert with detailed measurements of the flow-form interactions in a laboratory flume. The work presented herein focuses primarily on experiments conducted in a unique flow facility wherein high-resolution measurements of the turbulent flow field associated with the early stages of protodune development are obtained utilizing particle-image velocimetry (PIV) in a refractive-index-matched (RIM) environment. The RIM technique facilitates flow measurements extremely close to model surfaces as well as unimpeded optical access which are critical to understanding the flow-form coupling. A series idealized, fixed-bed models are fabricated to mimic the key morphological characteristics of early protodune development observed in the field, and the flow measurements associated with them are analyzed to reveal the mechanisms controlling the bedform dynamics.

**Tjebbe Hepkema**, Huib de Swart, and Henk Schuttelaars

**Dmitry Dukhovskoy**

Increasing Greenland discharge has contributed more than 5000 km^{3} of surplus fresh water to the Subpolar North Atlantic since the early 1990s. The volume of this freshwater anomaly is projected to cause freshening in the North Atlantic leading to changes in the intensity of deep convection and thermohaline circulation in the subpolar North Atlantic. This is roughly half of the freshwater volume of the Great Salinity Anomaly of the 1970s that caused notable freshening in the Subpolar North Atlantic. In analogy with the Great Salinity Anomaly, it has been proposed that, over the years, this additional Greenland freshwater discharge might have a great impact on convection driving thermohaline circulation in the North Atlantic with consequent impact on climate. Previous numerical studies demonstrate that roughly half of this Greenland freshwater anomaly accumulates in the Subpolar Gyre. However, time scales over which the Greenland freshwater anomaly can accumulate in the subpolar basins is not known. This study estimates the residence time of the Greenland freshwater anomaly in the Subpolar Gyre by approximating the process of the anomaly accumulation in the study domain with a first order autonomous dynamical system forced by the Greenland freshwater anomaly discharge. General solutions are obtained for two types of the forcing function. First, the Greenland freshwater anomaly discharge is a constant function imposed as a step function. Second, the surplus discharge is a linearly increasing function. The solutions are deduced by utilizing results from the numerical experiments that tracked spreading of the Greenland fresh water with a passive tracer. The residence time of the freshwater anomaly is estimated to be about 10–15 years. The main differences in the solutions is that under the linearly increasing discharge rate, the volume of the accumulated Greenland freshwater anomaly in the Subpolar Gyre does not reach a steady state. By contrast, solution for the constant discharge rate reaches a steady state quickly asymptoting the new steady state value for time exceeding the residence time. Estimated residence time is compared with the numerical experiments and observations.

**Tom Manzocchi**, Deirdre Walsh, Carneiro Marcus, Javier López-Cabrera, and Soni Kishan

Irrespective of the specific technique (variogram-based, object-based or training image-based) applied, geostatistical facies models usually use facies proportions as the constraining input parameter to be honoured in the output model. The three-dimensional interconnectivity of the facies bodies in these models increases as the facies proportion increases, and the universal percolation thresholds that define the onset of macroscopic connectivity in idealized statistical physics models define also the connectivity of these facies models. Put simply, the bodies are well connected when the model net:gross ratio exceeds about 30%, and because of the similar behaviour of different geostatistical approaches, some researchers have concluded that the same threshold applies to geological systems.

In this contribution we contend that connectivity in geological systems has more degrees of freedom than it does in conventional geostatistical facies models, and hence that geostatistical facies modelling should be constrained at input by a facies connectivity parameter as well as a facies proportion parameter. We have developed a method that decouples facies proportion from facies connectivity in the modelling process, and which allows systems to be generated in which both are defined independently at input. This so-called compression-based modelling approach applies the universal link between the connectivity and volume fraction in geostatistical modelling to first generate a model with the correct connectivity but incorrect volume fraction using a conventional geostatistical approach, and then applies a geometrical transform which scales the model to the correct facies proportions while maintaining the connectivity of the original model. The method is described and illustrated using examples representative of different geological systems. These include situations in which connectivity is both higher (e.g. fluid-driven injectite or karst networks) and lower (e.g. many depositional systems) than can be achieved in conventional geostatistical facies models.

**Peter Bossew**

The asymptotic shape of the marginal frequency distribution of geochemical variables has been proposed as indicator of multi-fractality. Transition into a certain statistical regime as inferred from the distribution function may be considered as criterion to delineate geochemical anomalies, including mineral resources or pollutants such as radioactive fallout or geogenic radon.

The argument is that asymptotic linearity in log-log scale, log(F(z)) = a - b log(z) as z→∞, b>0 a constant, indicates multi-fractality.

We discuss this with respect to two issues:

(1) What are the consequences of estimating the slope b for non-ergodic, in particular non-representative and preferential sampling schemes, as often the case in geochemical or pollution surveys?

(2) Frequently in geochemistry, multiplicative cascades are considered valid generators of multifractal fields, corroborated by observed f(α) functions and variograms (Matèrn or power, for low lags). This generator leads to marginally asymptotically (high cascade orders) log-normal distributions, which in log-log scale are asymptotically (high z) parabolic, not linear.

Theoretical aspects are addressed as well as examples given.

**Tommaso Alberti**, Giuseppe Consolini, Peter D. Ditlevsen, Reik V. Donner, and Virgilio Quattrociocchi

Several attempts have been made in characterizing the multiscale nature of fluctuations from nonlinear and nonstationary time series. Particularly, the study of their fractal structure has made use of different approaches like the structure function analysis, the evaluation of the generalized dimensions, and so on. Here we report on a different approach for characterizing phase-space trajectories by using the empirical modes derived via the Empirical Mode Decomposition (EMD) method. We show how the derived Intrinsic Mode Functions (IMFs) can be used as source of local (in terms of scales) information allowing us in deriving multiscale measures when looking at the behavior of the generalized fractal dimensions at different scales. This formalism is applied to three pedagogical examples like the Lorenz system, the Henon map, and the standard map. We also show that this formalism is readily applicable to characterize both the behavior of the Earth’s climate during the past 5 Ma and the dynamical properties of the near-Earth electromagnetic environment as monitored by the SYM-H index.

**Mark McGuinness**and Emma Greenbank

A Surtseyan volcanic eruption involves a bulk interaction between water and hot magma, mediated by the return of ejected ash. Surtsey Island, off the coast of Iceland, was born during such an eruption process in the 1940s. Mount Ruapehu in New Zealand also undergoes Surtseyan eruptions, due to its crater lake.

One feature of such eruptions is ejected lava bombs, trailing steam, with evidence that watery slurry was trapped inside them during the ejection process. Simple calculations indicate that the pressures developed due to boiling inside such a bomb should shatter it. Yet intact bombs are routinely discovered in debris piles. In an attempt to crack this problem, and provide a criterion for fragmentation of Surtseyan bombs, a transient mathematical model of the flashing of water to steam inside one of these hot erupted lava balls is developed, with a particular focus on the maximum pressure attained, and how it depends on magma and fluid properties. Numerical and asymptotic solutions provide some answers for volcanologists.

**Michal Kuraz**and Petr Mayer

Modeling the kinematic wave equation and sediment transport equation enables a deterministic approach for predicting surface runoff and resulting sediment transport. Both the kinematic wave equation and the sediment transport equation are first order differential equations. Moreover the kinematic wave equation is a quasilinear problem. In many engineering applications this set of equations is solved on one-dimensional representative cross-sections. However, a proper selection of representative cross-section(s) is cumbersome. On the other hand integrating this set of equations on real catchment topography yields difficulties for standard variational methods such as continous Galerkin method. These difficulties are two-fold (1) the nonlinearity of the kinematic wave, and (2) the absence of diffusion term, which acts as a stabilization term for convection-diffusion equation. In a theory, the Peclet number of numerical stability reaches infinity.

In this contribution we will focus on a stable numerical approximation of this convection-only problem using least square method. With this method we are able to reliably solve both the kinematic wave equation and the sediment transport equation on computational domains representing real catchment topography. Several examples representing real-world scenarios will be given.

**Elena Malinovskaya**and Otto Chkhetiani

One of the important characteristics of the wind process of dust removal is a critical or threshold wind velocity [1]. Saltating flow grows with increasing of the effective roughness [2] that affecting shear stress and friction velocity [3]. The drag coefficient increases depending on the density of the coating by particles of the surface [4]. The location of particles in the aeolian structure, their size and relative position determine their resistance to wind influence. Aeolian structures change the structure of flows and the balance of mass transfer of particles deposited and rising from the surface [5]. The surface microstructures and ripples significantly affect of sand removal.

The flow of particles with a size of 100 μm on the surface has been considered using the OPEN FOAM with LES model. The calculation area has sizes of 5x5x2 mm. For the velocity at the upper boundary, 2.8 m/s select in accordance with the experimental data [6]. It should be noted that with a relative increase in the distance between pairs of particles and a change in the level of the upper surface, the pressure difference between the base and top of the particle increases by 10-30 percents. Depending on the distance between the particles, the buoyant force acting from the side of the air flow, the critical velocity, and the departure velocity of the particle also change. When the distances between the surfaces of the particles are close to its size, the buoyant force is greater than the adhesion and gravity forces. As a result, areas with different probability for the sand removal by wind, due to which, in particular, the occurrence of aeolian ripples occurs.

The average critical velocity increases when moving up the windward slope of the dune [7, 8]. This phenomenon is possibly associated with the influence of ripples on the air flow. The flow around of the micro-ripples with a height of 0.1-1 mm was considered for air flow velocity of 2-4 m/s at a height of 1-2 cm. The addition of supplementary elements of heterogeneity at the apex near the rough surface of the streamlined aeolian structure leads to a displacement of the separation point of the ascending flows. Also we have a change in the length of the recirculation zone and the time intervals of the strengthening of the wind at the apex, which was observed in [6].

The study was supported by the RFBR project 19-05-50110 and partial support of the program of the Presidium of the Russian Academy of Sciences No. 12.

1. Shao Y. Physics and modeling of wind erosion. Springer.2008.p.452.

2. Martin R.L., Kok J.F. J.Geophys.Res.2018.123(7).1546-1565.

3. Turpin C et al. Earth Surf. Proc. and Land.2010.35(12). 1418-1429.

4. Yang X.I.A. et al. J. Fluid Mech.2019.880. 992-1019.

5. Luna M.C.M.M. et al. Geomorph.2011.129(3-4). 215-224.

6. Semenov O.E. Introduction to experimental meteorology and climatology of the sand storms. Almaty. 2011. p.580 (in Russian).

7. Neuman C.M.K et al. Sediment. 2000. 47(1). 211-226.

8. Malinovskaya E.A. Izv. Atmos. Oceanic Phys. 2019. 55(2). 86-92.

## NP2 – Dynamical Systems and Stochastic Approaches in Geosciences

### NP2.1 – Nonlinear Multiscale and Stochastic Dynamics of the Earth System

**Yuchen Ma**and William Peltier

Global coupled climate modeling requires the representation of multiple widely separated scales in each model component. In the ocean component, the separation of scales is especially dramatic as small scale turbulence exerts significant control on the global scale overturning circulation. The importance of this control is demonstrated in the context of analyses of the Dansgaard-Oeschger oscillation of Marine Isotope Stage 3 (MIS 3; see Peltier and Vettoretti, 2014)). In the University of Toronto version of CCSM4 water column diapycnal diffusivity is represented using the KPP parameterization of Large et al (1994). This includes explicit contributions due to double diffusion processes which demonstrably play an important role in determining the period of the D-O oscillation itself.

We have developed a DNS-based methodology to test the accuracy of the doubly diffusive contributions to KPP. High-resolution turbulence data sets have been produced based upon two different models: the “unbounded gradient model” and the “interface model” with depth-dependent temperature and salinity gradients. By fitting the vertical fluxes in the unbounded gradient model (for equilibrium states) as a function of density ratio (the governing non-dimensional parameter) we derive a functional form on the basis of which KPP can be revised. By applying the revised scheme to the interface model we demonstrate that the local fluxes predicted agree well with those from the numerical simulations. The difference between this new parametrization scheme and KPP demonstrates that KPP may significantly overestimate the diffusivities for both heat and salt at low-density ratio.

**Vinicius Beltram Tergolina**, Stefano Berti, and Gilmar Mompean

When studying the life cycle of phytoplankton frequently one is interested in the survival or death conditions of a population (bloom/no bloom). These dynamics have been studied extensively in the literature through a range of modelling scenarios but in summary the main factors affecting the vertical dynamics are: Water column mixing intensity, solar energy distribution, nutrients availability and predatory activity. The later two can be represented by different biological models whereas the vertical mixing is usually parameterized by a diffusive process. Even though turbulence has been recognized as a paramount factor in the survival dynamics of sinking phytoplankton species, dealing with the multi scale nature of turbulence is a formidable challenge from the modelling point of view. In addition, convective motions are being recognized to play a role in the survival of phytoplankton throughout winter stocking. With this in mind, in this work we revisit a theoretically appealing model for phytoplankton vertical dynamics with turbulent diffusivity and numerically study how large-scale fluid motions affect its survival and extinction conditions. To achieve this and to work with realistic parameter values, we adopt a kinematic flow field to account for the different spatial and temporal scales of turbulent motions. The dynamics of the population density are described by a reaction-advection-diffusion model with a growth term proportional to sun light availability. Light depletion is modelled accounting for water turbidity and plankton self-shading; advection is represented by a sinking speed and a two-dimensional, multiscale, chaotic flow. Preliminary results show that under appropriate conditions for the flow, our model reproduces past results based on turbulent diffusivity. Furthermore, the presence of large scale vortices (such as those one might expect during winter convection) seems to hinder survival, an effect that is partially mitigated by turbulent diffusion.

**Ana M. Mancho**, Guillermo Garcia-Sanchez, and José Antonio Jimenez-Madrid

The European Commission has invested in developing services such as the Copernicus Marine Environment Monitoring Services that offer opportunities to new downstream applications. This presentation describes the development of monitoring services in coastal areas at the submesoscale, by addressing synergies between different available marine technologies and products such as satellite images, autonomous surface and underwater vehicles, drone images, downscaled hydrodynamic models, etc, that get inspired in recent success cases [1, 2]. In particular ongoing efforts will be discussed that address the operational implementation of these tools for the management of marine pollution in harbors and coasts with a focus in the hydrodynamic modelling aspects.

Support is acknowledged from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 821922 (IMPRESSIVE) and from Fundacion Biodiversidad and European Commission (BEWATS).

**References**

[1] A. G. Ramos, V. J. García-Garrido, A. M. Mancho, S. Wiggins, J. Coca, S. Glenn, O. Schofield, J. Kohut, D. Aragon, J. Kerfoot, T. Haskins, T. Miles, C. Haldeman, N. Strandskov, B. Allsup, C. Jones, J. Shapiro. Lagrangian coherent structure assisted path planning for transoceanic autonomous underwater vehicle missions. Sci. Rep. 8, 4575 (2018).

[2] V. J. Garcia-Garrido, A. Ramos, A. M. Mancho, J. Coca, S. Wiggins. A dynamical systems perspective for a real-time response to a marine oil spill. Mar. Pollut. Bull. 112, 201-210, (2016).

**Yang Gao**, Francois G Schmitt, Jianyu Hu, and Yongxiang Huang

Turbulence or turbulence-like phenomena are ubiquitous in nature, often showing a power-law behavior of the Fourier power spectrum in either spatial or temporal domains. This power-law behavior is due to interactions among different scales of motion, and to the absence of characteristic scale among several scale ranges. It can be further interpreted in the framework of turbulent cascade with movements on continuous range of scales. The power-law feature and the associate cascade picture are vitally important to our understanding of the ocean and atmosphere dynamics. In this work, we consider the China France Oceanography SATellite (CFOSAT) data in the general framework of ocean and atmosphere multi-scale dynamics. We apply both Fourier power spectrum analysis and second-order structure-function analysis, used in the fields of turbulence, to extract multiscale information from the wind speed (WS) and significant wave-height (Hs) data provided by CFOSAT project. The data analyzed here are along track data spatially collected from 29^{th} July to 31^{th} December 2019. The measured Fourier power spectrums for both WS and Hs illustrate a dual power-law behavior respectively from 5 to 25 km, and 30 to 500 km with measured scaling exponents β close to 2 and 5/3. The measured second-order structure-functions confirm the existence of the dual power-law behavior. The corresponding measured scaling exponents ζ(2) close to 1 and 2/3 for the spatial scales mentioned above. Our preliminary results confirm the relevance of using multiscale statistical tools and turbulent theory to characterize the large-scale movements of both ocean and atmosphere.

**Maristella Berta**, Annalisa Griffa, Lorenzo Corgnati, Marcello Magaldi, Carlo Mantovani, Helga Huntley, Andrew Poje, and Tamay Ozgokmen

The dynamics of the near-surface ocean currents result from the nonlinear interaction of simultaneous mechanisms at different scales. Of these, the submesoscale range (a few hundred meters to 10 km, hours to a few days) remains particularly challenging to observe directly, due to the high variability in both time and space. In this study, the scale-dependence of kinematic properties (divergence, vorticity and strain) in the submesoscale range, as well as their response to atmospheric forcing, is investigated in two distinct geographic regions: the Ligurian (NW-Mediterranean) Sea and the Northern Gulf of Mexico. The two applications are characterized by different dynamics, and the estimates of kinematic properties are derived from distinctly different observational approaches: in situ GPS drifters observations and remote HF radar data.

The Ligurian Sea application is based on HF radar measurements obtained for the JERICO-NEXT (Joint European Research Infrastructure network for Coastal Observatory – Novel European eXpertise for coastal observaTories) and IMPACT (Port Impact on Protected Marine Areas: Cooperative Cross-Border Actions) projects. Surface current measurements span 40 km off the coast with 1.5 km resolution, available every hour. The velocity fields are used to estimate the kinematic properties with an Eulerian approach, which allows the identification of structures such as eddies and jets of the order of a few km. We focus in particular on the response of the submesoscales to an extreme wind event that was captured by the observations. The deformation of the spatial structures suggests nonlinear interactions with the wind forcing, and the kinematic properties' magnitudes are almost doubled (exceeding the Coriolis parameter, f).

In the Gulf of Mexico, velocity observations are available from a series of massive, nearly simultaneous drifter releases conducted by CARTHE (Consortium for Advanced Research of Transport of Hydrocarbons in the Environment). Drifter triplets are analysed to estimate the kinematic properties of the flow at scales between 100 m and 5 km over a time scale of a day. Results show that the mean kinematic properties increase in magnitude with decreasing scales, with winter values generally higher than summer ones. For winter flows, vorticity and divergence distributions have more substantial tails of values multiple times the Coriolis paramater f at scales up to 2 km, while in the summer such large values are restricted to smaller scales (100-500 m).

The Lagrangian estimates of kinematic properties obtained in the Gulf of Mexico were also compared to Eulerian estimates from concurrent X-band radar measurements, showing good correlation and validating the comparison across observational methods. Moreover, the scale-dependence of the kinematic properties from drifter triplets was found to be consistent with turbulence scaling laws evaluated as two-particle statistics. We conclude that the kinematic properties metric provides a robust complementary methodology to characterize submesoscales and can be used both with Lagrangian and Eulerian observations.

**Jih-Wang Aaron Wang**and Prashant Sardeshmukh

Despite decades of development, global atmospheric models continue to have trouble in capturing the -5/3 slope of the atmospheric mesoscale kinetic energy (KE) spectrum suggested by conventional turbulence theory and upper tropospheric aircraft observations in the 1980s (e.g., Nastrom and Gage 1986). We have approached this issue in two ways: 1) How certain can we be that the “real” spectrum has a -5/3 slope? and 2) Are turbulent cascades the only determinants of the mesoscale spectrum? To address the first issue, especially in light of the vastly greater number of upper-air observations that have been analyzed since the 1980s, we have examined the 200-hPa KE spectra in several high-resolution global reanalysis datasets, including the NCEP GFS (resolution T1534 and T254), ERA-Interim (T255), ERA5 (T639), and JRA-55 (T319) datasets. We find that the mesoscale portions of the global spectra are highly mutually inconsistent. This is primarily because the global mesoscale KE has a large contribution from the KE in convective regions, which differs greatly among the various reanalyses. There is thus indeed some ambiguity concerning the slope of the “true” mesoscale spectrum.

To address the second issue, especially given the sensitivity of the reanalysis spectra to representations of convection and damping in the reanalysis models, we assessed the sensitivity of the model spectra in two ways: (a) by stochastically perturbing the physical tendencies and (b) by decreasing the hyper-viscosity coefficient, in large ensembles of 10-day forecasts made with the NCEP GFS (T254) model. Both changes increased the mesoscale KE and decreased the steep spectral slope. The impact of the stochastic physics varied with the specified length scale of the stochastic perturbations.

Our conclusions about issues 1) and 2) raised above are that (1) we do not really know the “true” mesoscale KE spectrum, and (2) model KE spectra are sensitive to and can be manipulated by stochastically perturbing the parameterized physical tendencies and tuning the horizontal diffusion in a model. It may therefore be misleading for modelers to pursue the -5/3 slope of the Nastrom-Gage spectrum.

**Yayun Zheng**

An abrupt climatic transition could be triggered by a single extreme event, and an alpha-stable non-Gaussian Levy noise is regarded as a type of noise to generate such extreme events. In contrast with the classic Gaussian noise, a comprehensive approach of the most probable transition path for systems under alpha-stable Levy noise is still lacking. We develop here a probabilistic framework, based on the nonlocal Fokker-Planck equation, to investigate the maximum likelihood climate change for an energy balance system under the influence of greenhouse effect and Levy fluctuations. We find that a period of the cold climate state can be interrupted by a sharp shift to the warmer one due to larger noise jumps with low frequency. Additionally, the climate change for warming 1.5 degree under an enhanced greenhouse effect generates a step-like growth process. These results provide important insights into the underlying mechanisms of abrupt climate transitions triggered by a Levy process.

**Georg Gottwald**and Federica Gugole

We employ the framework of the Koopman operator and dynamic mode decomposition to devise a computationally cheap and easily implementable method to detect transient dynamics and regime changes in time series. We argue that typically transient dynamics experiences the full state space dimension with subsequent fast relaxation towards the attractor. In equilibrium, on the other hand, the dynamics evolves on a slower time scale on a lower dimensional attractor. The reconstruction error of a dynamic mode decomposition is used to monitor the inability of the time series to resolve the fast relaxation towards the attractor as well as the effective dimension of the dynamics. We illustrate our method by detecting transient dynamics in the Kuramoto-Sivashinsky equation. We further apply our method to atmospheric reanalysis data; our diagnostics detects the transition from a predominantly negative North Atlantic Oscillation (NAO) to a predominantly positive NAO around 1970, as well as the recently found regime change in the Southern Hemisphere atmospheric circulation around 1970.

**Vera Melinda Galfi**, Lesley de Cruz, Valerio Lucarini, and Sebastian Schubert

We analyze linear perturbations of a coupled quasi-geostrophic atmosphere-ocean model based on Covariant Lyapunov Vectors (CLVs). CLVs reveal the local geometrical structure of the attractor, and point into the direction of linear perturbations applied to the trajectory. Thus they represent a link between the geometry of the attractor and basic dynamical properties of the system, and they are physically meaningful. We compute the CLVs based on the so-called Ginelli method using the tangent linear version of the quasi-geostrophic atmosphere-ocean model MAOOAM (Modular Arbitrary-Order Ocean-Atmosphere Model). Based on the CLVs, we can quantify the contribution of each model variable on each scale to the development of linear instabilities. We also study the changes in the structure of the attractor - and, consequently, in the basic dynamical properties of our system - as an effect of the ocean-atmopshere coupling strength and the model resolution.

**Artur Prugger**and Jens Rademacher

Large scale motions in geophysical fluid models are often characterised by linear waves, which are obtained by linearising the equations. But there are also many explicit solutions of the fully nonlinear equations when posed the full space. The exact solutions we are investigating often characterise Rossby waves, since they are in geostrophic balance. They also can be compositions of waves, some are interacting with each other and some do not, showing wave interactions as explicit solutions in the fully nonlinear problem.

In this talk I will briefly introduce the idea behind these explicit nonlinear waves and show some of their properties, and their occurrence in different fluid models in extended domains.

As an application, we especially focus on a rotating shallow water model with simplified backscatter. In this case one finds not only geostrophic explicit solutions, but also ageostrophic ones. Moreover, here energy accumulates in selected scales due to the backscatter terms and causes exponentially and unboundedly growing ageostrophic nonlinear waves. This also relates to instability of coexisting stationary waves and is an instance of the role of nonlinear waves in energy transfer, and illustrates their role in preventing energy equidistribution for general data.

**James Annan**, Julia Hargreaves, Thorsten Mauritsen, and Bjorn Stevens

We examine what can be learnt about climate sensitivity from variability in the surface air temperature record over the instrumental period, from around 1880 to the present. While many previous studies have used the trend in the time series to constrain equilibrium climate sensitivity, it has recently been argued that temporal variability may also be a powerful constraint. We explore this question in the context of a simple widely used energy balance model of the climate system. We consider two recently-proposed summary measures of variability and also show how the full information content can be optimally used in this idealised scenario. We find that the constraint provided by variability is inherently skewed and its power is inversely related to the sensitivity itself, discriminating most strongly between low sensitivity values and weakening substantially for higher values. As a result of this, is only when the sensitivity is very low that the variability can provide a tight constraint. Our results support the analysis of variability as a potentially useful tool in helping to constrain equilibrium climate sensitivity, but suggest caution in the interpretation of precise results.

**Anna von der Heydt**and Peter Ashwin

The equilibrium climate sensitivity (ECS) is widely used as a measure for possible future global warming. It has been determined from a wide range of climate models, observations and palaeoclimate records, however, it still remains relatively unconstrained. In particular, large values of warming as a consequence of atmospheric greenhouse gas increase cannot be excluded, with some of the most recent state-of-the-art climate models (CMIP6) supporting (much) more warming than previous generations of climate models. Moreover, a number of tipping elements have been identified within the climate system, some of which may affect the global mean temperature. Therefore, it is interesting to explore how the climate systems response (e.g. ECS) behaves when the system is close to a tipping point.

A climate state close to a tipping point will have a degenerate linear response to perturbations, which can be associated with extreme values of the equilibrium climate sensitivity (ECS). In this talk we contrast linearized ('instantaneous') with fully nonlinear geometric ('two-point') notions of ECS, in both presence and absence of tipping points. For a stochastic energy balance model of the global mean surface temperature with two stable regimes, we confirm that tipping events cause the appearance of extremes in both notions of ECS. Moreover, multiple regimes with different mean sensitivities are visible in the two-point ECS. We confirm some of our findings in a physics-based multi-box model of the climate system.

**Reference**

P. Ashwin and A. S. von der Heydt (2019), Extreme Sensitivity and Climate Tipping Points, J. Stat. Phys. **370**, 1166–24. http://doi.org/10.1007/s10955-019-02425-x.

**Valerio Lembo**, Gabriele Messori, Rune Graversen, and Valerio Lucarini

The atmospheric meridional energy transport in the Northern Hemisphere midlatitudes is mainly accomplished by planetary and synoptic waves. A decomposition into wave components highlights the strong seasonal dependence of the transport, with both the total transport and the contributions from planetary and synoptic waves peaking in winter. In both winter and summer months, poleward transport extremes primarily result from a constructive interference between planetary and synoptic motions. The contribution of the mean meridional circulation is close to climatology. Equatorward transport extremes feature a mean meridional equatorward transport in winter, while the planetary and synoptic modes mostly transport energy poleward. In summer, a systematic destructive interference occurs, with planetary modes mostly transporting energy equatorward and synoptic modes again poleward. This underscores that baroclinic conversion dominates regardless of season in the synoptic wave modes, whereas the planetary waves can be either free or forced, depending on the season.

**Lenin Del Rio Amador**and Shaun Lovejoy

From hourly to decadal time scales, atmospheric fields are characterized by two scaling regimes: at high frequencies the weather, with fluctuations increasing with the time scale, and at low frequencies, macroweather with fluctuations decreasing with scale, the transition between the two at *τ _{w}*. This transition time is the lifetime of planetary structures and is therefore close to the deterministic predictability limit of conventional numerical weather prediction models. While it is thus the outer scale of deterministic weather models, conversely, it is the inner scale of stochastic macroweather models.

Here we explore the spatial dependence of this transition time. Starting at the surface (2m temperature) we found that the monthly average temperature falls in the macroweather regime for almost any location in the globe, except for parts of the tropical ocean where *τ _{w }*∼ 1 - 2 years. As we increase in altitude, the dependence of

*τ*with the location becomes more homogeneous and above 850mb

_{w}*τ*< 1 month almost everywhere. The longer tropical ocean transition scales are presumably the deterministic outer scales of the “ocean weather” regime.

_{w}Knowledge of *τ _{w}* is fundamental for stochastic macroweather forecasting. Such forecasting is based on symmetries, primarily the power-law behavior of the fluctuations that implies a huge memory that can be exploited for forecasts up to several years. In addition, there is another approximate symmetry called “statistical space-time factorization” relating spatial and temporal statistics. Finally, while weather regime temperature fluctuations are highly intermittent, in macroweather the intermittency is much lower, fluctuations are quasi Gaussian.

The Stochastic Seasonal and Interannual Prediction System (StocSIPS^{[1,2]}) is a stochastic data-driven model that exploits these symmetries to perform macroweather (long-term) forecasts. Compared to traditional global circulation models (GCM), it has the advantage of forcing predictions to converge to the real-world climate (not the model climate). It extracts the internal variability (weather noise) directly from past data and does not suffer from model drift. Some other practical advantages include much lower computational cost, no need for downscaling and no ad hoc postprocessing.

We show that StocSIPS can predict monthly average surface temperature (nearly) to its stochastic predictability limits. Using monthly to annual lead time hindcasts, we compare StocSIPS predictions with those from the CanSIPS^{[3]} GCM. Beyond a month, and especially over land, StocSIPS generally has higher skill. For regular StocSIPS forecasts, see http://www.physics.mcgill.ca/StocSIPS/.

**References**

^{[1]} Del Rio Amador, L. and Lovejoy, S. (2019) Clim Dyn, **53**: 4373. https://doi.org/10.1007/s00382-019-04791-4

^{[2]} Lovejoy, S., Del Rio Amador, L., Hébert, R. (2017) In Nonlinear Advances in Geosciences, A.A. Tsonis ed. Springer Nature, 305–355 DOI: 10.1007/978-3-319-58895-7

^{[3]} Merryfield WJ, Denis B, Fontecilla JS, Lee WS, Kharin S, Hodgson J, Archambault B (2011) Rep., 51pp, Environment Canada.

**Francois Schmitt**, Hussein Yahia, Joel Sudre, Véronique Garçon, and Guillaume Charria

Oceanic fields display a large variability over large temporal and spatial scales. One way to characterize such variability, borrowed from the field of turbulence, is to consider scaling regimes and multi-scaling properties.

He we use 2D power spectral analysis as well as 2D structure functions <X(M)-X(N)^{q}>=F(q,d(M,N)), between tow points M and N belonging to the region of interest. By performing statistics with respect to the distance d(M,N), one may extract the scaling property of the 2D field, for a range of distances L_{min}<d<L_{max}, of the form F(q,d)=d^{ζ(q)}. This approach can be used even for irregular images (having missing values due to cloud coverage) or for part of images in order to estimate the statistical heterogeneity of different zones of a given image.

In the framework of the French CNRS/IMECO project, we consider MODIS Aqua SST images, in France (English Channel versus Gascogne Golf) and in Chili (Eastern Boundary Upwelling System). We illustrate the use of the 2D structure function analysis for different part of these images and also different times. Scaling ranges and also scaling exponents are compared. To take into account the anisotropy of some of these zones, an anisotropic version of the 2D structure functions is also used.

**Arne Bendinger**, Johannes Karstensen, Julien Le Sommer, Aurélie Albert, and Fehmi Dilmahamod

Mesoscale eddies play an important role in lateral property fluxes. Observational studies often use sea level anomaly maps from satellite altimetry to estimate eddy statistics (incl. eddy kinetic energy). Recent findings suggest that altimetry derived eddy characteristics may suffer from the low spatial resolution of past and current satellite-tracks in high-latitude oceans associated with small Rossby radii. Here we present results of an eddy reconstruction based on a nonlinear damping Gauss-Newton optimisation algorithm using ship based current profiler observations from two research expeditions in the Labrador Sea in 2014 and 2016. Overall we detect 14 eddies with radii ranging from 7 to 35 km.

In order to verify the skill of the reconstruction we used the submesoscale permitting NATL60 model (1/60°) as a reference data set. Spectral analysis of the horizontal velocity implies that the mesoscale regime is well represented in NATL60 compared with the observations. The submesoscale regime in the model spectra shows deviations to the observations at scales smaller than 10km near the ocean surface. The representation of the submesoscale flow further decreases in the model with increasing depth.

By subsampling the NATL60 model velocities along artificial ship tracks, applying our eddy reconstruction algorithm, and comparing the results with the full model field, a skill assessment of the reconstruction is done. We show that the reconstruction of the eddy characteristics can be affected by the location of the ship track through the velocity field.

In comparison with the observed eddies the NATL60 eddies have smaller radii and higher azimuthal velocities and thus are more nonlinear. The inner core velocity structure for observations and NATL60 suggests solid body rotation for 2/3 of the radius. The maximum azimuthal velocity may deviate by up to 50% from solid body rotation.

The seasonality of the submesoscale regime can be seen in the data as the power spectrum is reduced from spring to summer in both the ship-based measurements and model.

**Guillaume Charria**, Sébastien Theetten, Adam Ayouche, Coline Poppeschi, Joël Sudre, Hussein Yahia, Véronique Garçon, and François Schmitt

The Bay of Biscay and the English Channel, in the North-eastern Atlantic, are considered as a natural laboratory to explore the coastal dynamics at different spatial and temporal scales. In those regions, the coastal circulation is constrained by a complex topography (e.g. varying width of the continental shelf, canyons), river runoffs, strong tides and a seasonally contrasted wind-driven circulation.

Based on different numerical model experiments (from 400m to 4km spatial resolution, from 40 to 100 sigma vertical layers using 3D primitive equation ocean models), different features of the Bay of Biscay and English Channel circulation are assessed and explored. Both spatial (submesoscale and mesoscale) and temporal (from hourly to monthly) scales are considered. Modelled spatial scales, with a specific focus on the variability of fine scale features (e.g. fronts, filaments, eddies), are compared with remotely sensed observations (i.e. Sea Surface Temperature). Different methodologies as singularity and Lyapunov exponents allow describing fine scales features and are applied on both modelled and observed datasets. For temporal scales, in situ high frequency surface temperature measurements from coastal moorings (from COAST-HF observing network) provide a reference for the temporal variability to be modelled. Exploring differences in the temporal scales (from an Empirical Mode Decomposition) advises on the efficiency of our coastal modelling approach.

This result overview in the Bay of Biscay and the English Channel aims illustrating the input of coastal modelling activities in understanding multi-scale interactions (spatial and temporal).

**Ashwita Chouksey**, Xavier Carton, and Jonathan Gula

In recent years, the oceanographic community has devoted considerable interest to the study of SCVs (Submesoscale Coherent Vortices, i.e. vortices with radii between 2-30 km, below the first internal radius of deformation); indeed, both mesoscale and submesoscale eddies contribute to the transport and mixing of water masses and of tracers (active and passive), affecting the heat transport, the ventilation pathways and thus having an impact on the large scale circulation.

In different areas of the ocean, SCVs have been detected, via satellite or in-situ measurements, at the surface or at depth. From these data, SCVs were found to be of different shapes and sizes depending on their place of origin and on their location. Here, we will concentrate rather on the SCVs at depth.

In this study, we use a high resolution simulation of the North Atlantic ocean with the ROMS-CROCO model. In this simulation, we also identify the SCVs at different depths and densities; we analyse their site and mechanism of generation, their drift, the physical processes conducting to this drift and their interactions with the surrounding flows. We also quantify their physical characteristics (radius, thickness, intensity/vorticity, bias in polarity: cyclones versus anticyclones). We provide averages for these characteristics and standard deviations.

We compare the model results with the observational data, in particular temperature and salinity profiles from Argo floats and velocity data from currentmeter recordings.

This study is a first step in the understanding of the formation, occurrences and structure of SCVs in the North Atlantic Ocean, of help to improve their in-situ sampling.

**Yongxiang Huang**, Yang Gao, Qianguo Xing, Francois Schmitt, and Jianyu Hu

Algal blooms, also known as ‘red tide’, are extremely harmful to the marine ecosystem since they infuse toxins into seawater and stifle oxygen in the water columns. Visually, they demonstrate rich patterns in spatial due to the interaction between the ocean current and the wind. Using the satelliate remote sensing data provided by the Chinese satellite Gaofeng 1, we first derive a normalized difference vegetation index (NDVI), which can be used to separate efficiently different types of cases, e.g., no algae bloom (NAB), macro algae bloom (MAB), and phytoplankton algae bloom (PAB), etc. The classical structure-function analysis is performed. Our preliminary results confirm the existence of the power-law behavior on the spatial scale range from 100 m to 400 m for the case of MAB. The corresponding scaling exponents are close to the ones of the classical passive scalar in three-dimension hydrodynamic turbulence. It suggests that the MAB could be treated as a passive scalar, which leads to not only a better understanding of the dynamics of algal blooms, but also a challenge of the modelling.

**Xia Zhang**, Liang Chen, Zhuguo Ma, and Yanhong Gao

The parameterization of surface exchange coefficients (C_{h}) representing land–atmosphere coupling strength plays a key role in land surface modeling. Previous studies have found that land–atmosphere coupling in land surface models (LSMs) is overestimated, which affects the predictability of weather and climate evolution. To improve the representation of land–atmosphere interactions in LSMs, this study investigated the dynamic canopy-height-dependent coupling strength in the offline Noah LSM with multiparameterization options (Noah-MP) when applied to China. Comparison with the default Noah-MP LSM showed the dynamic scheme significantly improved the C_{h} calculations and realistically reduced the biases of simulated surface energy and water components against observations. It is noteworthy that the improvements brought by the dynamic scheme differed across land cover types. The scheme was found superior in reproducing the observed C_{h} as well as surface energy and water variables for short vegetation (grass, crop, and shrub), while the improvement for tall canopy (forest) was found not significant, although the estimations were reasonable. The improved version benefits from the treatment of the roughness length for heat. Overall, the dynamic coupling scheme markedly affects the simulation of land–atmosphere interactions, and altering the dynamics of surface coupling has potential for improving the representation of land–atmosphere interactions and thus furthering LSM development.

**Christian Franzke**, Lichao Yang, and Zuntao Fu

Precipitation is an important meteorological variable which is critical for weather risk assessment. For instance, intense but short precipitation events can lead to flash floods and landslides. Most statistical modelling studies assume that the occurrence of precipitation events is based on a Poisson process with exponentially distributed waiting times while precipitation intensities are typically described by a gamma distribution or a mixture of two exponential distributions. Here, we show by using hourly precipitation data over the United States that the waiting time between precipitation events is non-exponentially distributed and best described by a fractional Poisson process. A systematic model selection procedure reveals that the hourly precipitation intensities are best represented by a two-distribution model for about 90% of all stations. The two-distribution model consists of (a) a generalized Pareto distribution (GPD) model for bulk precipitation event sizes and (b) a power-law distribution for large and extreme events. Finally, we analyse regional climate model output to evaluate how the climate models represent the high-frequency temporal structure of U.S. precipitation. Our results reveal that these regional climate models fail to accurately reproduce the power-law behaviour of intensities and severely underestimate the long durations between events.

**Abdel Hannachi**, Thomas Önskog, and Christian Franzke

The North Atlantic Oscillation (NAO) is the dominant mode of climate variability over the North Atlantic basin and has a significant impact on seasonal climate and surface weather conditions. This is the result of complex and nonlinear interactions between many spatio-temporal scales. Here, the authors study a number of linear and nonlinear models for a station-based time series of the daily winter NAO index. It is found that nonlinear autoregressive models including both short and long lags perform excellently in reproducing the characteristic statistical properties of the NAO, such as skewness and fat tails of the distribution and the different time scales of the two phases. As a spinoff of the modelling procedure, we are able to deduce that the interannual dependence of the NAO mostly affects the positive phase and that timescales of one to three weeks are more dominant for the negative phase. The statistical properties of the model makes it useful for the generation of realistic climate noise.

**Naiming Yuan**, Wenlu Wu, Fenghua Xie, and Yanjun Qi

Long-term persistence (LTP) and multifractality in river runoff fluctuations have been well recognized over the recent decades, but the origins of these characteristics are still under debate. In this study, runoff and precipitation data from China are analyzed using detrended fluctuation analysis (DFA) and its generalized version, multifractal detrended fluctuation analysis (MF-DFA). By comparing the results between runoff and the nearby precipitation data, we find the multifractal behaviors in river runoff may be propagated from the nearby precipitation data, but the LTP is not inherited from precipitation. The LTP in river runoff may arise from the spatial aggregation effect, as it is closely related with the catchment area, especially for stations with large catchment areas. These findings are based on data from China, which was not analyzed systematically due to the poor data availability. Since the existence of LTP and multifractality makes the runoff change not completely random, one should further introduce these characteristics into hydrological models, for improved water managements and better estimations of hazard risks.

**Carlos Pires**and Abdel Hannachi

The El-Niño index behaves as a nonlinear and non-Gaussian stochastic process. A well-known characteristic is its positive skewness coming from the occurrence of stronger episodes of El-Niño than of La Niña. Here, we use the period 1870-2018 of the standardized El-Niño index x(t), sampled in trimesters to analyze the spectral origin of the bicorrelation: sk(t1,t2)=E[x(t)x(t+t1)x(t+t2)] and skewness sk(0,0). For that, we estimate the two-dimensional Fourier transform of sk(t1,t2) or bispectrum B(f1,f2). Its sum over bi-frequencies (f1,f2) equals the skewness (0.45 in our case). Positive and negative bispectrum peaks are due to phase locking of frequency triplets: (f1,f2,f1+f2), contributing to extreme El-Niños and La Niñas respectively. Moreover, the most significant positive and/or negative bispectrum regions are rather well localized in the bispectrum domain. Here, we propose a partition of the El Niño signal into a set of band-pass spectrally separated components whose self and cross interactions can explain the broad structure of bispectrum. In the simplest case where the signal is decomposed into a fast and a slow component (with a cutoff frequency of (1/2.56) cycles/yr.), we verifty that slow-slow interactions (or phase locking) explain most of La-Niñas, particularly at the frequency triplet (1/4.9, 1/15 and 1/3.7 cycles/yr) whereas the fast-slow interactions explain most of El Niños, particularly at the frequency triplet (1/4.9, 1/4.9 and 1/2.5 cycles/yr). In order to simulate this stochastic behavior, we calibrate a set of nonlinearly coupled oscillators (Auto-regressive processes, forced by self and cross quadratic component terms), one for each component. In the case of weak cross-component interactions, and thus weak nonlinearity, the coupling coefficients between spectral-band components are proportional to the corresponding cross-skewnesses, which represent good first guesses in the calibration of the model parameters. The predictability of the model is then assessed, in particular for the anticipation of big El Niños and la Niñas. The authors would like to acknowledge the financial support FCT through project **UIDB/50019/2020 – IDL.**

**Shaobo Qiao**

Using observations and model simulations, the impact of the November Eurasian (EU) teleconnection on the following January Arctic Oscillation (Arctic Oscillation) and the possible mechanisms are investigated in this study. We found that the positive (negative) phase of the November EU pattern favors the negative (positive) phase of AO during the subsequent January, and both the stratosphere-troposphere interactionsand the tropospheric Blocking High (BH) activity anomalies over the Euro-Atlantic sector play an important role in their connections. When the EU pattern is positive (negative) phase in November, the increased (decreased) vertical wave activity over Eurasia and North Atlantic gradually weakens (enhances) the Stratospheric polar vortex (SPV)from November to the following early January, which is then conducive to a downward propagation of positive height anomalies from the stratosphere to troposphere. On the other hand, due to the persistent stronger (weaker) and southward (northward) shifted storm tracks over the Euro-Atlantic sector from November to the following early January, the BH activities over this region are significantly decreased (increased) during the same period, whichthen contributesto positive (negative) height anomalies over the Arctic via the propagation of a zonal wave number 1-3. As both the SPVand BHanomalies over the Euro-Atlantic sector reach the maximum around the late December-early January, the resultant equivalent barotropic AO dipole patterndevelops and finally establishes during the following January.These results are useful for the predictability of the middle winter climate

**Meng Zou**

Using hindcast and forecastdata from the National Centers for Environmental Prediction (NCEP) Climate Forecast System version 2 (CFSv2)for the period 1982-2017, we comprehensively assess the predictability of the climatology, interannual variability, and dominant modes of the wintertime 500 hPa geopotential height over Ural-Siberia (40-80°Nand 30-100°E). Although the climatic mean 500 hPa heightover Ural-Siberia simulated by NCEP CFSv2has a negative bias, especially over the eastern part of the region, NCEP CFSv2 well predicts the spatial distribution of the two major modes(EOF1 and EOF2) over this region 2 months in advance.The forecasting skill of the principal component (PC) of the two major modes,PC1 (PC2), is highest1 (0) month in advance, where the linear correlation coefficient between the predicted and observed time series reaches +0.36 (+0.67), exceeding the 95% confidence level. Conversely, the forecasting skill of PC1 (PC2) is very low 0 (1) month in advance. The main reason for the poorer(better) prediction of PC1 0 (1) month in advance is associated with a less (more) accurate response of the Eurasian teleconnection to SST anomalies over the southwestern Atlantic. For PC2, the better (poorer) prediction of PC2 0 (1) month in advance may be due to more (less) accurate responses of the stratospheric polar vortex and the Scandinavian teleconnection to the dipole SST anomalies over the North Pacific. These results are useful for evaluating the predictability of the East Asian winter climate.

**Maria Parfenova**and Igor I. Mokhov

The estimations of various factors influence on weather regimes formation in Russian regions in transitional (spring, fall) seasons are presented. Changes in those regimes comparing to the middle of 20th century are analyzed, considering atmospheric circulation features under the changes in meridional heat transfer and Rossby waves stationary modes. Using long-term observations of surface air temperature from several locations across Russia, the multimodal features of the probability density functions (PDF) in several decades of 20th and 21st centuries are identified. Focusing on surface temperature anomalies in transitional seasons, we examine the connection between the multimodal features of their PDFs and the nonlinear dependence of surface albedo on temperature during the formation and melting of snow cover. We investigate the impacts of other mechanisms that can facilitate these features, including blocking of zonal atmospheric transport in middle latitudes and formation of blocking anticyclones (blockings) and stationary Rossby waves.

**Benjamin Martinez-Lopez**

Sea surface temperature (SST) is the only oceanic parameter on which depend heat fluxes between ocean and atmosphere and, therefore, SST is one of the key factors that influence climate and its variability. Over the twentieth century, SSTs have significantly increased around the global ocean, warming that has been attributed to anthropogenic climate change, although it is not yet clear how much of it is related to natural causes and how much is due to human activities. A considerable part of available literature regarding climate change has been built based on the global or hemispheric analysis of surface temperature trends. There are, however, some key open questions that need to be answered and for this task estimates of long-term SST trend patterns represent a source of valuable information. Unfortunately, long-term SST trend patterns have large uncertainties and although SST constitutes one of the most-measured ocean variables of our historic records, their poor spatial and temporal sampling, as well as inhomogeneous measurements technics, hinder an accurate determination of long-term SST trends, which increases their uncertainty and, therefore, limit their physical interpretation as well as their use in the verification of climate simulations.

Most of the long-term SST trend patterns have been built using linear techniques, which are very usefull when they are used to extract information of measurements satisfying two key assumptions: linearity and stationarity. The global warming resulting of our economic activities, however, affect the state of the World Ocean and the atmosphere inducing changes in the climate that may result in oscillatory modes of variability of different frequencies, which may undergo non-stationary and non-linear evolutions. In this work, we construct long-term SST trend patterns by using non-linear techniques to extract non-linear, long-term trends in each grid-point of two available global SST datasets: the National Oceanic and Atmospheric Administration Extended Reconstructed SST (ERSST) and from the Hadley Centre sea ice and SST (HadISST). The used non-linear technique makes a good job even if the SST data are non-linear and non-stationary. Additionally, the nonlinearity of the extracted trends allows the use of the first and second derivative to get more information about the global, long-term evolution of the SST fields, favoring thus a deeper understanding and interpretation of the observed changes in SST. Particularly, our results clearly show, in both ERSST and HadISST datasets, the non-uniform warming observed in the tropical Pacific, which seems to be related to the enhanced vertical heat flux in the eastern equatorial Pacific and the strengthening of the warm pool in the western Pacific. By using the second derivative of the nonlinear SST trends, emerges an interesting pattern delimiting several zones in the Pacific Ocean which have been responded in a different way to the impose warming of the last century.

### NP2.3 – Extremes in Geophysical Sciences: Dynamics, Thermodynamics and Impacts

**Ángel G. Muñoz**

Cross-timescale interference involves linear and non-linear interactions between climate modes acting at multiple timescales (Muñoz et al., 2015, 2016, 2017; Robertson et al., 2015; Moron et al., 2015), and that are related to windows of opportunity for enhanced predictive skill (Mariotti et al., 2020), with relevant societal impacts (e.g., Doss-Gollin et al., 2018; Anderson et al., 2020). Using a simple mathematical model, reanalysis data and gridded observations, here we analyze plausible mechanisms for cross-timescale interference, describing conditions for coupling of oscillating modes and its impact on extreme rainfall occurrence and predictive skill. Concrete examples for Northeast North America and southern South America are discussed, as well as implications for climate model diagnostics.

**Giusy Fedele**, Thierry Penduff, Stefano Pierini, M. Carmen Alvarez-Castro, Alessio Bellucci, and Simona Masina

The Kuroshio Extension (KE) is the inertial meandering jet formed by the convergence of the Kuroshio and Oyashio currents in the Northern Pacific. It is widely mentioned in the literature that the KE variability is bimodal on the decadal time scale. The nature of this low frequency variability (LFV) is still under debate; some authors suggest that internal oceanic mechanisms play a fundamental role in the phenomenon but there is also evidence from the observations that the KE LFV is connected with changes in broader patterns of variability such as the Pacific Decadal Oscillation.

We first inspect the interplay between the ocean and the atmosphere in the KE by taking advantage of the OCCIPUT 1/4° model dataset: it consists in an ensemble of 50 global ocean–sea-ice hindcasts performed over the period 1960–2015 (hereafter OCCITENS), and in a one-member 330-yr climatological simulation (hereafter OCCICLIM). In this context, OCCITENS simulates both the intrinsic and forced variability and allows for their separation via ensemble statistics, while OCCICLIM simulates the "pure" intrinsic variability of the jet. We then explore some features of the KE dynamical system attractor in the quasi-autonomous (OCCICLIM) and nonautonomous (OCCITENS) regimes: we thus assess the KE predictability in the OCCIPUT dataset in order to better understand the ocean-atmosphere interactions and the source of the associated predictability.

Our analyses show that both oceanic and atmospheric drivers control the KE LFV variability. In this framework, the results suggest that the jet oscillates between the two intrinsic oceanic modes with transitions triggered by the atmosphere.

**Joakim Kjellsson**, Wonsun Park, Torge Martin, Eric Maisonnave, and Mojib Latif

We study how mesoscale air-sea interactions over the North Atlantic can influence weather extremes, e.g. heavy precipitation and wind storms, and the overall atmospheric circulation both locally and downstream in the midlatitudes. We use a global coupled climate model with a high-resolution North Atlantic grid (dx ~ 8 km) and an atmosphere model resolution of either 125 km or 25 km. The high-resolution North Atlantic grid allows the model to resolve the current systems and SST fronts associated with e.g. the Gulf Stream and North Atlantic Current. As air-sea fluxes of momentum, heat and freshwater are calculated on the atmosphere grid, spatial variations in fluxes associated with sharp SST fronts are much better represented when using the high-resolution atmosphere then when using the low-resolution model.

Preliminary results show that coupling to the high-resolution (dx ~ 25 km) rather than low-resolution (dx ~ 125 km) atmosphere model increases the intensity and variance of surface heat and freshwater fluxes over eddy-rich regions such as the Gulf Stream. As a result, the high-resolution model simulates more intense heavy precipitation events over most of the North Atlantic Ocean. We also show that more frequent coupling between the atmosphere and ocean components increases the intensity of the air-sea fluxes, in particular wind stress, which has a large impact on the ocean. More intense air-sea fluxes can provide more energy for cyclogenesis and we will discuss how the oceanic mesoscale, in particular in the eddy-rich regions, can alter the storm tracks and jet stream to influence extreme weather and the climate over Europe.

The coupled model comprises NEMO 3.6/LIM2 ocean and OpenIFS 40r1 atmosphere, and works by allowing the global OpenIFS model to send and receive fields from both a global coarse-resolution ocean grid and a refined grid over the North Atlantic grid via the OASIS3-MCT4 coupler. The ability to run these simulations is a very recent development and we will give a brief overview of the coupled modelling system and benefits of using regional grid refinement in coupled models.

**Kathrin Wehrli**, Mathias Hauser, and Sonia I. Seneviratne

The 2018 summer was unusually hot in large areas of the Northern Hemisphere and simultaneous heat waves on three continents led to major impacts to agriculture and society. The event was driven by the anomalous atmospheric circulation pattern during that summer and it was only possible in a climate with global warming. There are indications that in a future, warmer climate similar events might occur regularly, affecting major ‘breadbasket’ regions of the Northern Hemisphere.

This study aims to understand the role of climate change for driving the intensity of the 2018 summer and to explore the sensitivity to changing warming levels. Model simulations are performed using the Community Earth System Model to investigate storylines for the extreme 2018 summer given the observed atmospheric large-scale circulation but different levels of background global warming: no human imprint, the 2018 conditions, and different mean global warming levels (1.5°C, 2°C, 3°C, and 4°C). The storylines explore the consequences of the event in an alternative warmer or colder world and thus help to increase our understanding of the drivers involved. The results reveal a strong contribution by the present-day level of global warming and provide an outlook to similar events in a possible future climate.

**Antoine Blanc**, Juliette Blanchet, and Jean-Dominique Creutin

Large-scale circulations (LSCs) explain a significant part of Alpine precipitations. Characterizing circulations triggering heavy precipitation is usually done using weather-type classifications. A different characterization is implemented here, based on analogy using the atmospheric descriptors proposed in Blanchet et al 2018, 2019. These descriptors are both related to the dynamics of LSC and to their relative position in the atmospheric space. This work is applied to the Isère river catchment for the 1950-2011 period, considering a 3-day time step. The 500 hPa and 1000 hPa geopotential heights covering part of the western Europe are used separately to represent LSC. Two analogy criteria are investigated for constructing the atmospheric descriptors, namely TWS and RMSE.

Our results reveal that LSCs triggering heavy precipitation amounts correspond to strong geostrophic wind with quasi constant direction during the three days, corresponding to blocking situations in altitude. Moreover, those patterns of circulation are among the least singulars, and they show the highest degree of clustering in the atmospheric space. We interpret the latest results by the fact that heavy precipitation LSCs feature twin circulation patterns. In addition, the 500 hPa geopotential height appears to discriminate better heavy precipitation situations than the 1000 hPa one. Finally, our work points out the benefit of a combined use of TWS and RMSE. TWS gives information about the direction of geostrophic wind, while RMSE -combined with TWS- informs about its strength.

References:

Blanchet, J., Stalla, S., and Creutin, J.-D. (2018). Analogy of multi-day sequences of atmospheric circulation favoring large rainfall accumulation over the French Alps. Atmospheric Science Letters.

Blanchet, J., Creutin, J-D. (2019). Modelling rainfall accumulations over several days in the French Alps using low-dimensional atmospheric predictors based on analogy. Journal of Applied Meteorology and Climatology.

**Audrey Brouillet**and Sylvie Joussaume

Global warming is projected to intensify during the 21st century. This warming will be more readily perceived by human populations if it occurs rapidly and if it induces a thermal heat stress on the human body. Yet, only few studies investigate how climate change could be felt by future populations. Here we assess this possible perceived evolution between 1959 and 2100 only combining thermodynamic and statistical indicators. We analyse extremes of temperature (T_{99}) and simplified Wet-Bulb Globe Temperature (WBGT_{99}), a common heat stress index assessing the combined effect of elevated temperature and humidity on the human body. For each year of the period, we define the speed of change as a difference between two successive 20-year periods (i.e. with a moving baseline), and assess how these running changes emerge from each last 20-y inter-annual variability.

According to a subset of 12 CMIP5 Earth System Models and the RCP8.5 scenario, the change of T_{99} and WBGT_{99} will be twice as fast in the future compared to the current speed of change in the mid-latitudes, and by up to four times faster tropical regions such as Amazonia. Warming accelerations are thus similar for both T_{99} and WBGT_{99}. However, in tropical regions by 2080, the speed is projected to be 2.3 times larger than the recent inter-annual variability for WBGT_{99}, and only 1.5 to 1.8 times larger for T_{99}. Currently, speeds of change are only 0.2 to 0.8 times as large as the recent year-to-year variability for both metrics. We also show that 36% of the total world population will experience an emergent WBGT_{99} intensification in 2080, but only 15% of the population for T_{99}. According to future projections, the accelerated warming of future heat extremes will be more felt by populations than current changes, and this perceived change will be more severe for WBGT_{99} than for T_{99}, particularly in the tropics.

**Tim Cowan**, Gabi Hegerl, Luke Harrington, and Friederike Otto

The central United States experienced the hottest summers of the twentieth century in 1934 and 1936, with over 40 heat wave days and maximum temperatures surpassing 44°C at some locations like Kansas and Oklahoma. In fact, as of 2019, the summer of 1936 is still the hottest on record. The heat waves coincided with the decade-long Dust Bowl drought, that caused wide-spread crop failures, dust storms that penetrated to New York and considerable out-migration. In a very-large ensemble regional modelling framework, we show that greenhouse gas increases slightly enhanced the frequency and duration of the Dust Bowl heat waves, and would strongly enhance similar heat waves in the present day under current, further elevated greenhouse gas levels. Specifically, present-day atmospheric greenhouse gas forcing would reduce the return period of a rare (less than once in a century) heat wave summer as observed in 1936 to about 1-in 40-years, with further contribution by sea surface warming. Here, we show that a key driver of this elevated heat wave risk is the reduction in evaporative cooling and increase in sensible heating during dry springs and summers. Hence, we conclude that a warmer world is creating the potential for future extreme heat in moisture-limited regions, with potentially very damaging impacts.

**Theodore Shepherd**

Extreme climate events are invariably highly nonlinear, complex events, resulting from the confluence of multiple causal factors, and often quite singular. In any complex system there is a tension between analysis methods that respect the singularity of the extreme events at the price of statistical repeatability, and those that emphasize statistical repeatability at the price of nonlinearity and complexity; this dichotomy is found across all areas of science. In the climate context, the ‘storyline’ approach has emerged in recent years as a way of following the first of these two pathways. I will discuss how the storyline approach can be cast within the mathematical framework of causal networks, which provides a way to bridge between the storyline and probabilistic approaches. This also provides a way to interpret data in an appropriately conditional manner, thereby aiding model-measurement comparison.

**Matthias Röthlisberger**, Michael Sprenger, Emmanouil Flaounas, Urs Beyerle, and Heini Wernli

In the last decades, extremely hot summers (hereafter extreme summers) have challenged societies worldwide through their adverse ecological, economic and public health effects. In this study, extreme summers are identified at all grid points in the Northern Hemisphere in the upper tail of the July–August (JJA) seasonal mean 2-meter temperature (T2m) distribution, separately in ERA-Interim reanalyses and in 700 simulated years with the Community Earth System Model (CESM) large ensemble for present-day climate conditions. A novel approach is introduced to characterize the substructure of extreme summers, i.e., to elucidate whether an extreme summer is mainly the result of the warmest days being anomalously hot, or of the coldest days being anomalously mild, or of a general shift towards warmer temperatures on all days of the season. Such a statistical characterization can be obtained from considering so-called rank day anomalies for each extreme summer, that is, by sorting the 92 daily mean T2m values of an extreme summer and by calculating, for every rank, the deviation from the climatological mean rank value of T2m.

Applying this method in the entire Northern Hemisphere reveals spatially strongly varying extreme summer substructures, which agree remarkably well in the reanalysis and climate model data sets. For example, in eastern India the hottest 30 days of an extreme summer contribute more than 70% to the total extreme summer T2m anomaly, while the colder days are close to climatology. In the high Arctic, however, extreme summers occur when the coldest 30 days are substantially warmer than climatology. Furthermore, in roughly half of the Northern Hemisphere land area, the coldest third of summer days contribute more to extreme summers than the hottest third, which highlights that milder than normal coldest summer days are a key ingredient of many extreme summers. In certain regions, e.g., over western Europe and western Russia, the substructure of different extreme summers shows large variability and no common characteristic substructure emerges. Furthermore, we show that the typical extreme summer substructure in a certain region is directly related to the region’s overall T2m rank day variability pattern. This indicates that in regions where the warmest summer days vary particularly strongly from one year to the other, these warmest days are also particularly anomalous in extreme summers (and analogously for regions where variability is largest for the coldest days). Finally, for three selected regions, thermodynamic and dynamical causes of extreme summer substructures are briefly discussed, indicating that, for instance, the onset of monsoons, physical boundaries like the sea ice edge, or the frequency of occurrence of Rossby wave breaking, strongly determine the substructure of extreme summers in certain regions.

**Regina Rodrigues**, Andrea Taschetto, Alex Sen Gupta, and Gregory Foltz

In 2013/14 eastern South America experienced one of its worst droughts, leading to water shortages in São Paulo, the world’s fourth most populated city. This event was also responsible for a dengue fever outbreak that tripled the usual number of fatalities and reduced Brazilian coffee production leading to a global shortages and worldwide price increases. The drought was associated with an anomalous anticyclonic circulation off southeast South America that prevented synoptic systems reaching the region while inhibiting the development of the South Atlantic Convergence Zone and its associated rainfall. A concomitant and unprecedented marine heatwave also developed in the southwest Atlantic. Here we show from observations that such droughts and adjacent marine heatwaves have a common remote cause. Atmospheric blocking triggered by tropical convection in the Indian and Pacific oceans can cause persistent anticyclonic circulation that not only leads to severe drought but also generates marine heatwaves in the adjacent ocean. We show that increased shortwave radiation due to reduced cloud cover and reduced ocean heat loss from weaker winds are the main contributors to the establishment of marine heatwaves in the region. The proposed mechanism, which involves droughts, extreme air temperature over land and atmospheric blocking explains approximately 60% of the marine heatwave events in the western South Atlantic. We also identified an increase in frequency, duration, intensity and extension of marine heatwave events over the satellite period 1982–2016. Moreover, surface primary production was reduced during these events with implications for regional fisheries.

**Nathalie Schaller**, Jana Sillmann, Malte Müller, Reindert Haarsma, Wilco Hazeleger, Trine Jahr Hegdahl, Timo Kelder, Gijs van den Oord, Albrecht Weerts, and Kirien Whan

A physical climate storyline approach is applied to an autumn flood event caused by an atmospheric river in the West Coast of Norway. The aim is to demonstrate the value and challenges of higher spatial and temporal resolution in simulating impacts. The modelling chain used is the same as the one used operationally, to issue flood warnings for example. Its output is therefore familiar to many users, which we expect will facilitate stakeholder engagement. Two different versions of a hydrological model are run to show that on the one hand, the higher spatial resolution between the global and regional model is necessary to realistically simulate the high spatial variability of precipitation in such a mountainous region. On the other hand we also show that the intensity of the peak streamflow is only captured realistically with hourly data. The higher resolution regional atmospheric model is able to simulate the fact that with the passage of an atmospheric river, some valleys receive high amounts of precipitation and others not. However, the coarser resolution global model shows uniform precipitation in the whole region. Translating the event into the future leads to similar results: while in some catchments, a future flood might be 50% larger than a present one, in others no event occurs because the atmospheric river does not hit that catchment.

**Milan Palus**

The mathematical formulation of causality in measurable terms of predictability was given by the father of cybernetics N. Wiener [1] and formulated for time series by C.W.J. Granger [2]. The Granger causality is based on the evaluation of predictability in bivariate autoregressive models. This concept has been generalized for nonlinear systems using methods rooted in information theory [3,4]. The information-theoretic approach, defining causality as information transfer, has been successful in many applications and generalized to multivariate data and causal networks [e.g., 5]. This approach, rooted in the information theory of Shannon, usually ignores two important properties of complex systems, such as the Earth climate: the systems evolve on multiple time scales and their variables have heavy-tailed probability distributions. While the multiscale character of complex dynamics, such as air temperature variability, can be studied within the Shannonian framework [6, 7], the entropy concepts of Rényi and Tsallis have been proposed to cope with variables with heavy-tailed probability distributions. We will discuss how such non-Shannonian entropy concepts can be applied in inference of causality in systems with heavy-tailed probability distributions and extreme events, using examples from the climate system.

This study was supported by the Czech Science Foundation, project GA19-16066S.

[1] N. Wiener, in: E. F. Beckenbach (Editor), Modern Mathematics for Engineers (McGraw-Hill, New York, 1956)

[2] C.W.J. Granger, Econometrica 37 (1969) 424

[3] K. Hlaváčková-Schindler et al., Phys. Rep. 441 (2007) 1

[4] M. Paluš, M. Vejmelka, Phys. Rev. E 75 (2007) 056211

[5] J. Runge et al., Nature Communications 6 (2015) 8502

[6] M. Paluš, Phys. Rev. Lett. 112 (2014) 078702

[7] N. Jajcay, J. Hlinka, S. Kravtsov, A. A. Tsonis, M. Paluš, Geophys. Res. Lett. 43(2) (2016) 902–909

**Zhenghui Lu**, Naiming Yuan, Zhuguo Ma, Qing Yang, and Juergen Kurths

The different phases of the Pacific Decadal Oscillation (PDO) are a primary source of internal decadal climate variability which have distinct impacts on global climate and human society. However, obtaining a reliable prediction of the PDO phase transition is still challenging. Here, we employed the new technique of climate network analysis to uncover early warning signals prior to a PDO phase transition. An examination of cooperative behaviors in the PDO region revealed an enhanced signal that propagated from the western Pacific, passed through the Kuroshio extension (KE) and the subtropical oceanic frontal (STF) regions, and finally reached the northwest coast of the Americas. This signal captured all six of the PDO phase transitions from the 1890s to 2000s, with a warning time of 6.5±2.3 years in advance. It also underpinned the possible PDO phase transition at years around 2015, which may be triggered by the strong El Niño in 2014-2016.

**Frank Kwasniok**

Traditional extreme value analysis based on the generalised ex-

treme value (GEV) or generalised Pareto distribution (GPD) suffers

from two drawbacks: (i) Both methods are wasteful of data as only

block maxima or exceedances over a high threshold are taken into ac-

count and the bulk of the data is disregarded. (ii) Moreover, in the

GPD approach, there is no systematic way to determine the threshold

parameter. Here, all the data are fitted simultaneously using a gener-

alised exponential family model for the bulk and a GPD model for the

tail. At the threshold, the two distributions are linked together with

appropriate matching conditions. The model parameters are estimated

from the likelihood function of all the data. Also the threshold param-

eter can be determined via maximum likelihood in an outer loop. The

method is exemplified on wind speed data from an atmospheric model.

**Gabriele Messori**, Erica Madonna, Orly Lachmy, and Davide Farranda

Atmospheric jet streams are typically separated into primarily "eddy-driven", or "polar-front" jets and primarily "thermally-driven", or "subtropical" jets. Some regions also display “merged” jets, resulting from the (quasi) co-location of the regions of eddy generation with the subtropical jet. The different location and driving mechanisms of the two jet structures, plus the intermediate “merged” jet, issue from very different underlying mechanisms, and result in very different jet characteristics. Here, we link our understanding of the dynamical jet maintenance mechanisms, mostly issuing from conceptual or idealised models, to the phenomena observed in reanalysis data. We specifically focus on developing a unitary analysis framework, grounded in dynamical systems theory, which may be applied to both the model and reanalysis data and allow for direct intercomparison. Our results provide a proof-of-concept for using dynamical systems indicators to diagnose jet regimes in a versatile, conceptually intuitive and computationally efficient fashion.

**Davide Faranda**, and Dim Coumou

The Mediterranean (MED) basin is a hot-spot for climate change impacts. We present recently developed techniques derived from Dynamical System Theory to investigate long-term changes in compound hot-wet extremes over the MED. We use three reanalysis products, spanning a 40-year period from 1979 to 2018: ERA5, ERA-Interim and ERA5 10-member ensemble. From these datasets, we extract daily maximum temperature (degC) and total precipitation (mm), which we then use in the dynamical systems analysis.

Results show that the strength of the dynamical coupling between hot and wet extremes increased significantly at both annual and summer (June-August) timescales during the reanalysis period. This means that, regardless of changes in the occurrence of individual hot or wet extremes, joint occurrences may be becoming more frequent.

Compound hot-wet extremes mostly occur during July and August, and correspond to a low-pressure core over the Aegean Sea and the eastern MED. The increasing trends in compound extremes may be associated with surface MED warming. Such enhanced warming can therefore drive compound hot-wet extremes especially during the summer, when localised convection or mesoscale systems such as medicanes are responsible for extreme precipitation events.

**Flavio Pons**, Paolo De Luca, Gabriele Messori, and Davide Faranda

We propose a novel approach to the study of compound extremes, grounded in dynamical systems theory. Specifically, we present the co-recurrence ratio (α), which elucidates the dependence structure between maps by quantifying their joint recurrences. This approach is applied to daily climate extremes, derived from the ERA-Interim reanalysis over the 1979-2018 period. The analysis focuses on concurrent (i.e. same-day) wet (total precipitation) and windy (10m wind gusts) extremes in Europe and concurrent cold (2m temperature) extremes in Eastern North America and wet extremes in Europe. Results for wet and windy extremes in Europe, which we use as a test-bed for our methodology, show that α peaks during boreal winter. High αvalues correspond to wet and windy extremes in north-western Europe, and to large-scale conditions resembling the positive phase of the North Atlantic Oscillation (NAO). This confirms earlier findings which link the positive NAO to a heightened frequency of extra-tropical cyclones impacting north-western Europe, resulting in frequent wet and windy extremes. For the Eastern North America-Europe case, α extremes once again reflect concurrent climate extremes -- in this case cold extremes over North America and wet extremes over Europe. Our analysis provides detailed spatial information on regional hotspots for these compound extreme occurrences, and encapsulates information on their spatial footprint which is typically not included in a conventional co-occurrence analysis. We conclude that α successfully characterises compound extremes by reflecting the evolution of the associated meteorological maps. This approach is entirely general, and may be applied to different types of compound extremes and geographical regions.

**Jana Sillmann**, and Marion Haugen

Climate models aim to project future changes in important drivers of climate including atmosphere, oceans and ice, and their interactions. A comprehensive evaluation of climate models thus requires evaluation methods, or performance measures, that are flexible, specific and can address also extreme events. Climate models have traditionally been assessed by comparing summary statistics or point estimates that derive from the simulated model output to corresponding observed quantities using e.g. RMSE. However, it has been argued persuasively that probability distributions of model output need to be compared to the corresponding empirical distributions of observations or observation-based data products. Observation-based gridded datasets for climate extremes, despite having limitations, are particularly useful and necessary to assess model performance with respect to extremes. We discuss proper performance measures for comparing distributions of model output against corresponding distributions from data products that are flexible and robust enough to handle the particular aspects of extremes such as limited data availability. The new measures are applied to evaluate CMIP5 and CMIP6 projections of extreme temperature indices over Europe and North-America against the HadEX2 data set as well as the ERA5 and ERA-Interim reanalyses. Several models perform well to the extent that when compared to the HadEX2 data product, these models' performance is competitive with the performance of the reanalysis. While the model rankings vary with region, season and index, the model evaluation is robust against changes in the grid resolution considered in the analysis.

**Theophile caby**, Davide Faranda, Sandro Vaienti, and Pascal Yiou

We study the properties of recurrence of a smooth observable computed along a trajectory of a chaotic system near a particular value of interest . Using the framework of Extreme Value Theory, we are able to derive a limit law which is a Gumbel distribution whose parameters relate to the dimensions of the image measure. We show that this approach allows to have access to the fine structure of the attractor, by using directly observational data. In particular, we are able to compute local dimensions associated to the underlying attractor whenever the dimensionality of the observable is larger than the dimension of the attractor.

**Arielle Catalano**, Paul Loikith, and J. David Neelin

Under global warming, changes in extreme temperatures will manifest in more complex ways in locations where temperature distribution tails deviate from Gaussian. For example, uniform warming applied to a temperature distribution with a shorter-than-Gaussian warm tail would lead to greater exceedances in warm-side temperature extremes compared with a Gaussian distribution. Confidence in projections of future temperature extremes and associated impacts under global warming therefore relies on the ability of global climate models (GCMs) to realistically simulate observed temperature distribution tail behavior. This presentation examines the ability of the latest state-of-the-art ensemble of GCMs from the Coupled Model Intercomparison Project phase six (CMIP6) to capture historical global surface temperature distribution tail shape in hemispheric winter and summer seasons. Comparisons between the multi-model ensemble mean and a reanalysis product reveal strong agreement on coherent spatial patterns of longer- and shorter-than-Gaussian tails for the cold and warm sides of the temperature distribution, suggesting that CMIP6 is broadly capturing tail behavior for plausible physical and dynamical reasons. Most individual GCMs are also reasonably skilled at capturing historical tail shape on a global scale, but a division of the domain into sub-regions reveals considerable model and spatial variability. To explore potential mechanisms driving these differences, a back trajectory analysis examining patterns in the origin of air masses on days experiencing extreme temperatures is also discussed.

**Elsa Salzedas**

“Water in the solar system” EGU 2020

** **

**Abstract **

Elsa Salzedas, Afonso de Albuquerque High School, Guarda, Portugal. (elsasalzedas@aeaag.pt) or ( elsa.salzedas@gmail.com)

**The thermal waters of the Portugal Star Geopark - methods for understanding its origin and sustainable exploitation.**

- How do the hot springs surface?

- What is the relationship between hot springs and tectonic accidents?

- How long does it take a drop of thermal water to reach the surface after it has infiltrated the soil?

- What precautions should be taken for sustainable exploitation.

- What precautions should be taken for a contaminant free operation.

- How to make people aware of the importance of thermal waters and their sustainable exploitation.

The Afonso de Albuquerque school, my school, is in the area of Geopark Star, which is an excellent laboratory for understanding these phenomena thermal waters.

High school students show us these discoveries through very interesting schemes, photographs and models, which will be presented at the next GIFT EGU in Mai 2020 in Vienna, Austria.

**Peter Pfleiderer**, Aglaé Jézéquel, Juliette Legrand, Natacha Legrix, Jason Markantonis, and Edoardo Vignotto

In 2016, northern France experienced an unprecedented wheat crop loss. This extreme event was likely due to particular meteorological conditions, i.e. too few cold days in late autumn and an abnormally high precipitation during the spring season. The cause of this event is not fully understood yet and none of the most used crop forecast models were able to predict the event (Ben-Ari et al, 2018).

This work is motivated by two main questions: were the 2016 meteorological conditions the most extreme we could imagine under current climate? and what would be the worst case scenario we could expect that could lead to even worst crop loss? To answer these questions, instead of relying on computationally intensive climate model simulations, we use an analogue based importance sampling algorithm that was recently introduced into this field of research (Yiou and Jézéquel, 2019). This algorithm is a modification of a stochastic weather generator (SWG), which gives more weight to trajectories with more extreme meteorological conditions (here temperature and precipitation). This data driven technique constructs artificial weather events by combining daily observations in a dynamically realistic manner and in a relatively fast way.

This is the first application of SWGs to simulate warm winters and wet springs. We show that with some adjustments both (new) weather events can be adequately simulated with SWGs, highlighting the wide applicability of the method.

While the number of cold days in late autumn 2015 was close to the plausible maximum, our simulations of extreme spring precipitation show that considerably wetter springs than what was observed in 2016 are possible. Although the crop loss of 2016 is not fully understood yet, these results indicate that similar events with higher impacts could be possible.

**Jongsoo Shin**and Soon-Il An

Heatwaves are likely to occur more frequent, longer, and stronger due to the rise in CO2 concentrations. It is related to the change in the mean of a climate distribution, as well as through the change in variance. Mega-heatwaves, in particular, have a crucial impact on human health. Many studies are trying to understand the mechanisms of mega-heatwaves and also their characteristics included amplitude, duration, frequency. In spite of these efforts, researches are limited because of the small number of mega-heatwaves. In order to overcome these limitations, Earth system model should be needed. This study aims to figure out the comprehensive characteristics of mega-heatwaves using Community Earth System Model (CESM). First, the possibility of the occurrence of mega-heatwaves in preindustrial period in Europe was analyzed. Second, the relation between decadal climate variabilities and mega-heatwaves was investigated. In addition, changes in characteristics of mega-heatwaves were compared between preindustrial and present-day simulations.

**Wenguang Wei**, Zhongwei Yan, and Zhen Li

On the decadal time scales, while the influence of Pacific Decadal Oscillation (PDO) on total or average precipitation had been extensively studied, works about its influence on precipitation extremes were limited, especially lack of a global picture. Using two independent methods, nonstationary generalized extreme value (GEV) model which directly incorporates PDO index into its location parameter and moving GEV model which fits the annual extremes with a sliding window of 30 years and regresses the resulted changing location parameter onto the PDO index, we show that precipitation extremes over a large portion of stations are significantly affected by the PDO with stations in the Pacific Rim demonstrating distinct regional patterns. Over eastern China, the famous ‘southern flood and norther drought’ pattern corresponding to a positive PDO phase extends to extreme rainfalls; over Australia, a tri-polar pattern was revealed, in which the extremes over central Australia positively correlate with the PDO index and those over eastern and western Australia show a negative correlation; and the North America also demonstrates a dipole pattern, by which the northwest (southeast) experiences less (more) intense extreme rainfall in a PDO positive phase. Moreover, the western Europe and the large area between the Ural mountain and eastern Europe were discovered to hold a positive correlation with the PDO in their precipitation extremes. A comparative analysis to the local circulation controlling the precipitation extremes under different PDO phases further confirms the discovered relationships above. These findings have important implication for the future projection of extreme precipitation over related regions because the internal climate variability should be appropriately accounted for beyond the effects induced by global warming.

**Ji-Seon Oh**, Maeng-Ki Kim, Dae-Geun Yu, and Jeong Sang

In this study, we defined diagnostic indices to evaluate the CMIP6 models based on the heatwaves mechanisms of Korea presented in previous studies. Based on this, the simulation performance of the model was quantitatively evaluated using Relative Error (RE), Inter-annual Variability Skill-score (IVS), and Correlation Coefficient (CC). The REs in diagnostic indices are still large, especially in heat wave circulation index (HWCI) and Indian summer monsoon rainfall index (IMRI), which is mainly due to weak convective activity bias over the northwestern Pacific Ocean and the northwestern India. However, the IVSs in diagnostic indices have been improved overall in the CMIP6 compared to the CMIP5, especially in the IMRI. The CC between the daily maximum temperature (TMAX) and the diagnostic factors in the model is very higher in HWCI than other indices, indicating that the convective activity over the northwestern Pacific is very important in heat wave in Korea. As a result, the total ranking of the model performance for heatwaves in Korea suggested that EC-Earth3-Veg, EC-Earth3, and UKESM-1-0-LL ranked high in CMIP6.

This work was funded by the Korea Meteorological Administration Research and Development Program under Grant KMI(KMI2018-03410)

**M. Carmen Alvarez-Castro**, Silvio Gualdi, Pascal Yiou, Mathieu Vrac, Robert Vautard, Leone Cavicchia, David Gallego, Pedro Ribera, Cristina Pena-Ortiz, and Davide Faranda

Windstorms, extreme precipitations and instant floods seems to strike the Mediterranean area with increasing frequency. These events occur simultaneously during intense tropical-like Mediterranean cyclones. These intense Mediterranean cyclones are frequently associated with wind, heavy precipitation and changes in temperature, generating high risk situations such as flash floods and large-scale floods with significant impacts on human life and built environment. Although the dynamics of these phenomena is well understood, little is know about their climatology. It is therefore very difficult to make statements about the frequency of occurrence and its response to climate change. Thus, intense Mediterranean cyclones have many different physical aspects that can not be captured by a simple standard approach.

The first challenge of this work is to provide an extended catalogue and climatology of these phenomena by reconstructing a database of intense Mediterranean cyclones dating back up to 1969 using the satellite, the literature and reanalyses. Applying a method based on dynamical systems theory we analyse and attribute their future changes under different anthropogenic forcings by using future simulations within CMIP framework. Preliminary results show a decrease of the large-scale circulation patterns favoring intense Mediterranean cyclones in all the seasons except summer.

**Felix Soltau**, Matthias Hirt, Jessica Kelln, Sara Santamaria-Aguilar, Sönke Dangendorf, and Jürgen Jensen

In the past decades, severe so called ‚compound events‘ led to critical high water levels at the coasts of southern Africa and as a consequence to property damage and loss of human life. The co-occurrence of storm surges, wind waves, heavy precipitation and resulting runoff increases the risk of coastal flooding and exacerbates the impacts along the vulnerable southern African coasts (e.g. Couasnon et al. 2019). To mitigate these high-impacts, it is essential to understand the underlying processes and driving factors (Wahl et al. 2015). As compound flooding events at southern African coasts are dominated by wind waves, it is of great importance to investigate the regional wave climate to understand the wave forcing as well as the origin of the wave energy.

Wind waves around southern African coasts are affected by the complex interactions between the Agulhas current and the atmosphere. In the research project CASISAC* we analyse the present evolution of the Agulhas Current system and quantify its impact on the future regional wave climate. Ocean waves contributing to high sea levels can be generated offshore resulting in swell or closer to the coasts by strong onshore winds. To identify responsible atmospheric pressure fields that force high wind wave events we apply a hybrid approach: (1) linking south hemispheric pressure fields with offshore wave data using an artificial neural network and (2) determine the prevailing nearshore wave conditions by regional numerical wave propagation models (SWAN). By validating the modelled nearshore wave data from hindcast runs with wave buoy records, this approach allows us to predict future extreme wind wave events and thus potential flooding. In a next step, extreme value analysis is used to determine future return periods of extreme wave events. These results can contribute to the development of more reliable adaptive protection strategies for southern African coast.

*CASISAC (Changes in the Agulhas System and its Impact on Southern African Coasts: Sea level and coastal extremes) is funded by the German Federal Ministry of Education and Research (BMBF) through the project management of Projektträger Jülich PTJ under the grant number 03F0796C

Couasnon, Eilander, Muis, Veldkamp, Haigh, Wahl, Winsemius, Ward (2019): Measuring compound flood potential from river discharge and storm surge extremes at the global scale and its implications for flood hazard. In: Natural Hazards and Earth System Sciences, Discussion Paper, S. 1–24. DOI: 10.5194/nhess-2019-205, in review.

Wahl, Jain, Bender, Meyers, Luther (2015): Increasing risk of compound flooding from storm surge and rainfall for major US cities. In: Nature Climate Change 5 (12), S. 1093–1097. DOI: 10.1038/nclimate2736.

## NP4 – Time Series and Big data methods

### NP4.1 – Complex geoscientific time series: linear, nonlinear, and computer science perspectives

**Ray Huffaker**and Rafael Munoz-Carpena

The complex soil biome is a center piece in providing essential ecosystem services that humans rely on (carbon sequestration, food security, one-health interactions). Agricultural engineers and soil scientists are developing wireless sensor networks (WSN) that collect large/big data on the soil key state variables (water content, temperature, chemistry) to better understand the soil biome primary environmental drivers. The profession extracts information from WSN records with methods including soil-process modeling and artificial-intelligence (AI) algorithms. However, these approaches carry their own limitations. A recent review article faulted current soil-process modeling for inadequately detecting and resolving model structural (abstraction) errors. AI experts themselves caution against indiscriminant use of AI methods because of: a) problems including replication of past results due to inconsistent experimental methods; b) difficulty in explaining how a particular method arrives at its conclusions (the black box problem) and thus in correcting algorithms that learn ‘bad lessons’; and c) lack of rigorous criteria for selecting AI architectures. An alternative approach to address these limitations is to investigate new strategies for reducing large/big data problems into smaller, more interpretable causal abstractions of the soil system.

We develop an innovative data diagnostics framework—based on empirical nonlinear dynamics techniques from physics—that addresses the above concerns over soil-process modeling and AI algorithms. We diagnose whether WSN and other similar environmental large/big data are likely generated by dimension-reducing (i.e., dissipative) nonlinear dynamics. An n-dimensional nonlinear dynamic system is dissipative if long-term dynamics are bounded within m<<n dimensions, so that the problem of modeling long-term dynamics shrinks by the n-m inactive degrees of freedom. If so, long-term system dynamics can be investigated with relatively few degrees of freedom that capture the complexity of the overall system generating observed data. To make this diagnosis, we first apply signal processing to isolate structured variation (signal) from unstructured variation (noise) in large/big data time series records, and test signals for nonlinear stationarity. We resolve the structure of isolated signals by distinguishing between stochastic-forcing and deterministic nonlinear dynamics; reconstruct phase space dynamics most likely generating signals, and test the statistical significance of reconstructed dynamics with surrogate data. If the reconstructed phase space is dimension-reducing, we can formulate low-dimensional (phenomenological) ODE models to investigate nonlinear causal interactions between key soil environmental driving factors. When we do not diagnose dimension-reducing nonlinear real-world dynamics, then underlying dynamics are most likely high dimensional and the information-extraction problem cannot be shrunk without losing essential dynamic information. In this case, other high-dimensional analysis techniques like AI offer a better modeling alternative for mapping out interactions. Our framework supplies a decision-support tool for data practitioners toward the most informative and parsimonious information-extraction method—a win-win result.

We will share preliminary results applying this empirical framework to three soil moisture sensor time series records analyzed with machine learning methods in Bean, Huffaker, and Migliaccio (2018).

**Isabel Serra**, David Moriña, Pere Puig, and Álvaro Corral

Intense geomagnetic storms can cause severe damage to electrical systems and communications. this work proposes a counting process with Weibull inter-occurrence times in order to estimate the probability of extreme geomagnetic events. It is found that the scale parameter of the inter-occurrence time distribution grows exponentially with the absolute value of the intensity threshold defining the storm, whereas the shape parameter keeps rather constant. The model is able to forecast the probability of occurrence of an event for a given intensity threshold; in particular, the probability of occurrence on the next decade of an extreme event of a magnitude comparable or larger than the well-known Carrington event of 1859 is explored, and estimated to be between 0.46% and 1.88% (with a 95% confidence), a much lower value than those reported in the existing literature.

**Zahra Zali**, Frank Scherbaum, Matthias Ohrnberger, and Fabrice Cotton

Volcanic tremor is one of the most important signal in volcano seismology because of its potential to be a tool for forecasting eruptions and better understanding of underlying volcanic process. Despite different suggested mechanisms for volcanic tremor generation, the exact process of that is not well understood yet. This signal usually comes along with large number of earthquakes happening during unrest period that affect the shape and amplitude of tremor. A delicate signal processing is required to separate earthquakes and other transient signals from seismic waveform to derive a time series of volcanic tremor which can provide a new insight into tremor source investigations. Exploiting the idea of harmonic and percussive separation in musical signal processing we have developed a method to extract volcanic tremor and transient events from the seismic signal. By using the concept of periodicity as underlying generation process of tremor, we are able to extract the volcanic tremor signal based on the self similarity properties of spectra in time-frequency domain. The separation process results in two spectrograms representing repeating (long-lasting) and non-repeating (short-lived) patterns.

From the spectrogram of the repeating pattern we reconstruct the signal in time domain by adding the original spectrogram’s phase information, thus creating an modified version of the long-lasting tremor signal.

Further, we can derive a characteristic function for transient events by integrating the amplitude of the non-repeating spectrogram in each time frame. This function has non zero value in transient event instances and zero value in time periods devoid of such events. Considering transient events as earthquakes we apply an onset detector to time first arrivals of the transient signal by using the slope of the function. First we determine local maxima of the function showing good correspondence to even the tiniest transient signals. From the peak locations we calculate the slope of each point within a period of 6 seconds preceding each peak. The uncertainty of positive P peaks is up to 0.32 seconds which is equal to the hope size of the calculated spectrogram. The advantage of timing earthquakes through this method is the ability of detecting very low seismic events, although due to the small window size of short time Fourier transform the process is time consuming. The result of this study is promising, while further testing is on-going to validate the method as well as determine applications and limitations.

**Alexandre Hippert-Ferrer**, Yajing Yan, and Philippe Bolon

Time series analysis constitutes a thriving subject in satellite image derived displacement measurement, especially since the launching of Sentinel satellites which provide free and systematic satellite image acquisitions with extended spatial coverage and reduced revisiting time. Large volumes of satellite images are available for monitoring numerous targets at the Earth’s surface, which allows for significant improvements of the displacement measurement precision by means of advanced multi-temporal methods. However, satellite image derived displacement time series can suffer from missing data, which is mainly due to technical limitations of the ground displacement computation methods (e.g. offset tracking) and surface property changes from one acquisition to another. Missing data can hinder the full exploitation of the displacement time series, which can potentially weaken both knowledge and interpretation of the physical phenomenon under observation. Therefore, an efficient missing data imputation approach seems of particular importance for data completeness. In this work, an iterative method, namely extended Expectation Maximization - Empirical Orthogonal Functions (EM-EOF) is proposed to retrieve missing values in satellite image derived displacement time series. The method uses both spatial and temporal correlations in the displacement time series for reconstruction. For this purpose, the spatio-temporal covariance of the time series is iteratively estimated and decomposed into different EOF modes by solving the eigenvalue problem in an EM-like scheme. To determine the optimal number of EOFs modes, two robust metrics, the cross validation error and a confidence index obtained from eigenvalue uncertainty, are defined. The former metric is also used as a convergence criterion of the iterative update of the missing values. Synthetic simulations are first performed in order to demonstrate the ability of missing data imputation of the extended EM-EOF method in cases of complex displacement, gaps and noise behaviors. Then, the method is applied to time series of offset tracking displacement measurement of Sentinel-2 images acquired between January 2017 and September 2019 over Fox Glacier in the Southern Alps of New Zealand. Promising results confirm the efficiency of the extended EM-EOF method in missing data imputation of satellite image derived displacement time series.

**Ganesh Ghimire**, Navid Jadidoleslam, Witold Krajewski, and Anastasios Tsonis

Streamflow is a dynamical process that integrates water movement in space and time within basin boundaries. The authors characterize the dynamics associated with streamflow time series data from about seventy-one U.S. Geological Survey (USGS) stream-gauge stations in the state of Iowa. They employ a novel approach called visibility graph (VG). It uses the concept of mapping time series into complex networks to investigate the time evolutionary behavior of dynamical system. The authors focus on a simple variant of VG algorithm called horizontal visibility graph (HVG). The tracking of dynamics and hence, the predictability of streamflow processes, are carried out by extracting two key pieces of information called characteristic exponent, λ of degree distribution and global clustering coefficient, GC pertaining to HVG derived network. The authors use these two measures to identify whether streamflow process has its origin in random or chaotic processes. They show that the characterization of streamflow dynamics is sensitive to data attributes. Through a systematic and comprehensive analysis, the authors illustrate that streamflow dynamics characterization is sensitive to the normalization, and the time-scale of streamflow time-series. At daily scale, streamflow at all stations used in the analysis, reveals randomness with strong spatial scale (basin size) dependence. This has implications for predictability of streamflow and floods. The authors demonstrate that dynamics transition through potentially chaotic to randomly correlated process as the averaging time-scale increases. Finally, the temporal trends of λ and GC are statistically significant at about 40% of the total number of stations analyzed. Attributing this trend to factors such as changing climate or land use requires further research.

**Laurane Charrier**, Yajing Yan, Elise Koeniguer, Emmanuel Trouvé, Romain Millan, Jérémie Mouginot, and Anna Derkacheva

Glacier response to climate change results in natural hazards, sea level rise and changes in freshwater resources. To evaluate this response, glacier surface flow velocity constitutes a crucial parameter to study. Nowadays, more and more velocity maps at regional or global scales issued from satellite SAR and/or optical images tend to be available online or on-demand. Such amount of data requires appropriate data fusion strategies in order to generate displacement time series with improved precision and spatio-temporal coverage. The improved displacement time series can then be used by advanced multi-temporal analysis approaches for further physical interpretations of the phenomenon under observation. In this work, time series of Sentinel-2 (10~m resolution, every 5 days), Landsat-8 (15~m resolution, every 16 days) and Venus (5~m resolution, every 2 days) images acquired between January 2017 and September 2018, over the Fox glacier in the Southern Alps of New Zealand are investigated. Velocities are generated with an offset tracking technique using an automatic processing chain for every possible repeat cycles (2 days-100 days and 300 days to 400 days). Thousands of velocity maps are available, and they are subject to both uncertainty and data gaps. In order to produce a displacement time series as precise/complete as possible , we propose three fusion strategies: 1) use all the available Sentinel-2 displacement maps with different time spans. The goal is to construct a time series of displacement with respect to a common master by means of an inversion 2) take only Sentinel-2 displacement maps with as small time spans as possible, at the same time, keep as much as possible redundancy in the network to be able to construct a common master displacement time series by inversion 3) follow the previous strategy but use all available displacement maps from 3 sensors, with different temporal sampling and measurement precision taken into account. Afterwards, the common master displacement time series will be analysed by a data mining approach in order to extract unusual spatio-temporal patterns in the time series.

**Yunzhong Shen**, Fengwei Wang, and Qiujie Chen

Since a time series is usually incomplete, the missing data are usually interpolated before employing singular spectrum analysis (SSA). We develop a new SSA for processing incomplete time series based on the property that an original time series can be reproduced from its principal components which are then estimated based on minimum norm criterion. When an incomplete time series is polluted by multiplicative noise, we first convert the multiplicative noise to additive noise by multiplying the signal estimate of the time series, then process the time series with weighted SSA, where the weight factor is determined according to the variance of additive noise, since the converted additive noise is heterogeneous. The proposed SSA approach is employed to process the real incomplete time series data of suspended-sediment concentration from San Francisco Bay compared to the traditional SSA and homomorphic log-transformation SSA approach. The first 10 principal components derived by our proposed SSA approach can capture more of the total variance and with less fitting error than traditional SSA approach and homomorphic log-transformation SSA approach. Furthermore, the results from the simulation cases conform that our proposed SSA outperform both traditional and homomorphic log-transformation SSA approaches.

**Adrian Odenweller**and Reik Donner

The quantification of synchronization phenomena of extreme events has recently aroused a great deal of interest in various disciplines. Climatological studies therefore commonly draw on spatially embedded climate networks in conjunction with nonlinear time series analysis. Among the multitude of similarity measures available to construct climate networks, Event Synchronization and Event Coincidence Analysis (ECA) stand out as two conceptually and computationally simple nonlinear methods. While ES defines synchrony in a data adaptive local way that does not distinguish between different time scales, ECA requires the selection of a specific time scale for synchrony detection.

Herein, we provide evidence that, due to its parameter-free structure, ES has structural difficulties to disentangle synchrony from serial dependency, whereas ECA is less prone to such biases. We use coupled autoregressive processes to numerically study the sensitivity of results from both methods to changes of coupling and autoregressive parameters. This reveals that ES has difficulties to detect synchronies if events tend to occur temporally clustered, which can be expected from climate time series with extreme events exceeding certain percentiles.

These conceptual concerns are not only reproducible in numerical simulations, but also have implications for real world data. We construct a climate network from satellite-based precipitation data of the Tropical Rainfall Measuring Mission (TRMM) for the Indian Summer Monsoon, thereby reproducing results of previously published studies. We demonstrate that there is an undesirable link between the fraction of events on subsequent days and the degree density at each grid point of the climate network. This indicates that the explanatory power of ES climate networks might be hampered since trivial local properties of the underlying time series significantly predetermine the final network structure, which holds especially true for areas that had previously been reported as important for governing monsoon dynamics at large spatial scales. In contrast, ECA does not appear to be as vulnerable to these biases and additionally allows to trace the spatiotemporal propagation of synchrony in climate networks.

Our analysis rests on corrected versions of both methods that alleviate different normalization problems of the original definitions, which is especially important for short time series. Our finding suggest that careful event detection and diligent preprocessing is recommended when applying ES, while this is less crucial for ECA. Results obtained from ES climate networks therefore need to be interpreted with caution.

**Reik Donner**and Jaqueline Lekscha

Analysing palaeoclimate proxy time series using windowed recurrence network analysis (wRNA) has been shown to provide valuable information on past climate variability. In turn, it has also been found that the robustness of the obtained results differs among proxies from different palaeoclimate archives. To systematically test the suitability of wRNA for studying different types of palaeoclimate proxy time series, we use the framework of forward proxy modelling. For this, we create artificial input time series with different properties and compare the areawise significant anomalies detected using wRNA of the input and the model output time series. Also, taking into account results for general filtering of different time series, we find that the variability of the network transitivity is altered for stochastic input time series while being rather robust for deterministic input. In terms of significant anomalies of the network transitivity, we observe that these anomalies may be missed by proxies from tree and lake archives after the non-linear filtering by the corresponding proxy system models. For proxies from speleothems, we additionally observe falsely identified significant anomalies that are not present in the input time series. Finally, for proxies from ice cores, the wRNA results show the best correspondence with those for the input data. Our results contribute to improve the interpretation of windowed recurrence network analysis results obtained from real-world palaeoclimate time series.

**Holger Lange**, Michael Hauhs, Katharina Funk, Sebastian Sippel, and Henning Meesenburg

We analyze time series from several forested headwater catchments located adjacent to each other in the Bramke valley, Harz mountains (Germany) which are monitored since decades for hydrology, hydrochemistry and forest growth. The data sets include meteorological variables, runoff rates, streamwater chemical concentrations, and others. The basic temporal resolution is daily for hydrometeorology and two-weekly for streamwater chemistry (in addition, standing biomass of a Norway spruce stand is measured every couple of years).

A model was calibrated and run for the streamflow from one of the catchments, based on precipitation, temperature and (simulated) evapotranspiration of the growing trees, to elucidate the effect of forest growth on catchment hydrology.

The catchments exhibit long-term changes and spatial gradients related to atmospheric deposition, management and changing climate. After providing a short multivariate summary of the dataset, we present several nonlinear metrics suitable to detect and quantify subtle changes and to describe different behavior, both between different variables from the same catchment, as well as for the same variable across catchments. The methods include, but are not limited to: Tarnopolski analysis, permutation entropy and complexity, q- and α-complexities, and Horizontal Visibility Graphs.

The detection of these changes is remarkable, because linear trends have already been removed prior to analysis. Hence, their presence reflects intrinsic changes in the patterns of the time series. The metrics also allow for a detailed model evaluation from a nonlinear perspective.

An important methodological aspect is the temporal resolution of the time series. We investigate the scaling behavior of the nonlinear metrics through aggregation or decimation to coarser resolutions and conclude on what the scaling behavior may imply for inverse (hydrological) modelling tasks.

**Stephane Vannitsem**, and Hugues Goosse

Dynamical dependence between key observables and the surface mass balance (SMB) over Antarctica is analyzed in two historical runs performed with the MPI‐ESM‐P and the CESM1‐CAM5 climate models. The approach used is a novel method allowing for evaluating the rate of information transfer between observables that goes beyond the classical correlation analysis and allows for directional characterization of dependence. It reveals that a large proportion of significant correlations do not lead to dependence. In addition, three coherent results concerning the dependence of SMB emerge from the analysis of both models: (i) The SMB over the Antarctic Plateau is mostly influenced by the surface temperature and sea ice concentration and not by large‐scale circulation changes; (ii) the SMB of the Weddell Sea and the Dronning Maud Land coasts are not influenced significantly by the surface temperature; and (iii) the Weddell Sea coast is not significantly influenced by the sea ice concentration.

**Andreas Gerhardus**and Jakob Runge

Scientific inquiry seeks to understand natural phenomena by understanding their underlying processes, i.e., by identifying cause and effect. In addition to mere scientific curiosity, an understanding of cause and effect relationships is necessary to predict the effect of changing dynamical regimes and for the attribution of extreme events to potential causes. It is thus an important question to ask how, in cases where controlled experiments are not feasible, causation can still be inferred from the statistical dependencies in observed time series.

A central obstacle for such an inference is the potential existence of unobserved causally relevant variables. Arguably, this is more likely to be the case than not, for example unmeasured deep oceanic variables in atmospheric processes. Unobserved variables can act as confounders (meaning they are a common cause of two or more observed variables) and thus introduce spurious, i.e., non-causal dependencies. Despite these complications, the last three decades have seen the development of so-called causal discovery algorithms (an example being FCI by Spirtes et al., 1999) that are often able to identify spurious associations and to distinguish them from genuine causation. This opens the possibility for a data-driven approach to infer cause and effect relationships among climate variables, thereby contributing to a better understanding of Earth's complex climate system.

These methods are, however, not yet well adapted to some specific challenges that climate time series often come with, e.g. strong autocorrelation, time lags and nonlinearities. To close this methodological gap, we generalize the ideas of the recent PCMCI causal discovery algorithm (Runge et al., 2019) to time series where unobserved causally relevant variables may exist (in contrast, PCMCI made the assumption of no confounding). Further, we present preliminary applications to modes of climate variability.

**Mirko Stumpo**, Giuseppe Consolini, Tommaso Alberti, and Virgilio Quattrociocchi

The fundamental question what causes what has always been the motivating motto for natural sciences, being the study of causality a crucial point for characterizing dynamical relationships. In the framework of complex dynamical systems, both linear statistical tools and Granger causality models drastically fail to detect causal relationships between time series, while a powerful model-free statistical framework is offered by the information theory.

Here we discuss how to deal with the problem of measuring causal information in non-stationary complex systems by considering a local estimation of the information-theoretic functionals via an ensemble-based statistics. Then, its application for investigating the dynamical coupling and relationships between the solar wind and the Earth’s magnetosphere is also presented.

**Antonio Cicone**, Angela Stallone, Massimo Materassi, and Haomin Zhou

Nonlinear and nonstationary signals are ubiquitous in real life. Their time–frequency analysis and features extraction can help in solving open problems in many fields of research. Two decades ago, the Empirical Mode Decomposition (EMD) algorithm was introduced to tackle highly nonlinear and nonstationary signals. It consists of a local and adaptive data–driven method which relaxes several limitations of the standard Fourier transform and the wavelet Transform techniques, yielding an accurate time-frequency representation of a signal. Over the years, several variants of the EMD algorithm have been proposed to improve the original technique, such as the Ensemble Empirical Mode Decomposition (EEMD) and the Iterative Filtering (IF).

The versatility of these techniques has opened the door to their application in many applied fields, like geophysics, physics, medicine, and finance. Although the EMD– and IF–based techniques are more suitable than traditional methods for the analysis of nonlinear and nonstationary data, they could easily be misused if their known limitations, together with the assumptions they rely on, are not carefully considered. Here we call attention to some of the pitfalls encountered when implementing these techniques. Specifically, there are three critical factors that are often neglected: boundary effects; presence of spikes in the original signal; signals containing a high degree of stochasticity. We show how an inappropriate implementation of the EMD and IF methods could return an artefact–prone decomposition of the original signal. We conclude with best practice guidelines for researchers who intend to use these techniques for their signal analysis.

**Mikhail Kanevski**, Federico Amato, and Fabian Guignard

The research deals with an application of advanced exploratory tools to study hourly spatio-temporal air pollution data collected by NABEL monitoring network in Switzerland. Data analyzed consist of several pollutants, mainly NO2, O3, PM2.5, measured during last two years at 16 stations distributed over the country. The data are considered in two different ways: 1) as multivariate time series measured at the same station (different pollutants and environmental variables, like temperature), 2) as a spatially distributed time series of the same pollutant. In the first case, it is interesting to study both univariate and multivariate time series and their complexity. In the second case, similarity between time series distributed in space can signify the similar underlying phenomena and environmental conditions giving rise to the pollution. An important aspect of the data is that they are collected at the places of different land use classes – urban, suburban, rural etc., which helps in understanding and interpretation of the results.

Nowadays, unsupervised learning algorithms are widely applied in intelligent exploratory data analysis. Well known tasks of unsupervised learning include manifold learning, dimensionality reduction and clustering. In the present research, intrinsic and fractal dimensions, measures characterizing the similarity and redundancy in data and machine learning clustering algorithms were adapted and applied. The results obtained give a new and important information on the air pollution spatio-temporal patterns. The following results, between others, can be mentioned: 1) some measures of similarity (e.g., complexity-independent distance) are efficient in discriminating between time series; 2) intrinsic dimension, characterizing the ensemble of monitoring data, is pollutant dependent; 3) clustering of time series observed can be interpreted using the available information on land use.

**Carmelo Juez**and Estela Nadal-Romero

Sediment transport is the major driver of changes in most catchments systems. Beyond landscape evolution and river geomorphology, sediment dynamics are an important component of a number of physical, chemical and biological processes in river basins. Sediments thus impact the ecology of rivers, sustainability of human infrastructure and basin level fluxes of nutrients and carbon. For this reason, it is important to understand the temporal sediment response of mountain catchments regarding precipitation and run-off. This response is not unique and features intra-annual, annual and multi-year scales components. In this research, we analyse a humid mountain badland area located in the Central Spanish Pyrenees. This typology of badlands is characterized by its non-linearity and non-stationary precipitation and run-off cycles, which ultimately lead to complex sediment dynamics and yields. Based on spectral frequency analysis and wavelet decomposition we were able to determine the dominant time scales of the local hydrological and sediment dynamics. Intra-annual and annual time scales were linked with the climatological characteristics of the catchment site. The multi-year response in the sediment yields reveals the importance of the sediment storage/depletion cycle of the catchment. The frequency and amplitude of precipitation, run-off and sediment yields fluctuations were accurately predicted with the spectral frequency analysis and wavelet decomposition technique used.

**Anastasia Nekrasova**and Vladimir Kossobokov

The observed variability of seismic dynamics of the Kamchatka Region is characterized in terms of several moving averages, including (i) seismic rate, (ii) the Benioff strain release, (iii) inter-event time, τ, and (iv) the USLE control parameter, η (where USLE stands for Unified Scaling Law for Earthquakes, i.e. a generalization of the Gutenberg-Richter relationship accounting for naturally fractal distribution of earthquake loci, which states that the distribution of inter-event times τ depends only on the value of variable η).

The variability of seismic dynamics have been evaluated and compared at each of four out of ten separate seismic focal zones of the Kamchatka region and the adjacent areas defined by Levina et al. (2013), i.e., (1) seismic focal zone of the Kuril and South Kamchatka, (2) the northern part of the Kamchatka seismic focal zone, (3) commander segment of the Aleutian arc; and (4) the continental region of Kamchatka. In particular, we considered all magnitude 3.5 or larger earthquakes in 1996-2019 available from open data catalog of the Kamchatka Branch of GS RAS, Earthquakes Catalogue for Kamchatka and the Commander Islands (1962–present) http://sdis.emsd.ru/info/earthquakes/catalogue.ph).

## NP5 – Predictability

### NP5.1 – Data Assimilation, Predictability, Error Identification and Uncertainty Quantification in Geosciences

**Luca Cantarello**, Onno Bokhove, Gordon Inverarity, Stefano Migliorini, and Steve Tobias

Operational data assimilation (DA) schemes rely significantly on satellite observations with much research aimed at their optimisation, leading to a great deal of progress. Here, we investigate the impact of the spatial-temporal variability of satellite observations for DA: is there a case for concentrating effort into the assimilation of small-scale convective features over the large-scale dynamics, or vice versa?

We conduct our study in an isentropic one-and-a-half layer model that mimics convection and precipitation, a revised and more realistic version of the idealised model based on the shallow water equations in [1,2]. Forecast-assimilation experiments are performed by means of a twin-setting configuration, in which pseudo-observations from a high-resolution nature run are combined with lower-resolution forecasts. The DA algorithm used is the deterministic Ensemble Kalman Filter (see [3]). We focus our research on polar-orbit satellites regarding emitted microwave radiation.

We have developed a new observation operator and a representative observing system in which both ground and satellite observations can be assimilated. The convection thresholds in the model are used as a proxy for cloud formation, clouds, and precipitation. To imitate the use of weighting functions in real satellite applications, radiance values are computed as a weighted sum with contributions from both layers. In the presence of clouds and/or precipitation, we model the response of passive microwave radiation to either precipitating or non-precipitating clouds. The horizontal resolution of satellite observations can be varied to investigate the impact of scale-dependency on the analysis.

New, preliminary results from experiments including both transverse jets and rotation in a periodic domain will be reported and discussed.

References:

[1] Kent, T., Bokhove, O., & Tobias, S. (2017). Dynamics of an idealized fluid model for investigating convective-scale data assimilation. Tellus A: Dynamic Meteorology and Oceanography, 69(1), 1369332.

[2] Kent, T. (2016). An idealised fluid model for convective-scale NWP: dynamics and data assimilation (Doctoral dissertation, PhD Thesis, University of Leeds).

[3] Sakov, P., & Oke, P. R. (2008). A deterministic formulation of the ensemble Kalman filter: an alternative to ensemble square root filters. Tellus A: Dynamic Meteorology and Oceanography, 60(2), 361-371.

**Jean-Michel Brankart**

Many practical applications involve the resolution of large-size inverse problems, without providing more than a moderate-size sample to describe the prior probability distribution. In this situation, additional information must be supplied to augment the effective dimension of the available sample, for instance using a covariance localization approach. In this study, it is suggested that covariance localization can be efficiently applied to an approximate variant of the Metropolis/Hastings algorithm, by modulating the ensemble members by the large-scale patterns of other members. Modulation is used to design a (global) proposal probability distribution (i) that can be sampled at a very low cost, (ii) that automatically accounts for a localized prior covariance, and (iii) that leads to an efficient sampler for the augmented prior probability distribution or for the posterior probability distribution. The resulting algorithm is applied to an academic example, illustrating (i) the effectiveness of covariance localization, (ii) the ability of the method to deal with nonlocal/nonlinear observation operators and non-Gaussian observation errors, (iii) the reliability, resolution and optimality of the updated ensemble, using probabilistic scores appropriate to a non-Gaussian posterior distribution, and (iv) the scalability of the algorithm as a function of the size of the problem. The codes are openly available from github.com/brankart/ensdam.

**Milija Zupanski**

High-dimensional ensemble data assimilation applications require error covariance localization in order to address the problem of insufficient degrees of freedom, typically accomplished using the observation-space covariance localization. However, this creates a challenge for vertically integrated observations, such as satellite radiances, aerosol optical depth, etc., since the exact observation location in vertical does not exist. For nonlinear problems, there is an implied inconsistency in iterative minimization due to using observation-space localization which effectively prevents finding the optimal global minimizing solution. Using state-space localization, however, in principal resolves both issues associated with observation space localization.

In this work we present a new nonlinear ensemble data assimilation method that employs covariance localization in state space and finds an optimal analysis solution. The new method resembles “modified ensembles” in the sense that ensemble size is increased in the analysis, but it differs in methodology used to create ensemble modifications, calculate the analysis error covariance, and define the initial ensemble perturbations for data assimilation cycling. From a practical point of view, the new method is considerably more efficient and potentially applicable to realistic high-dimensional data assimilation problems. A distinct characteristic of the new algorithm is that the localized error covariance and minimization are global, i.e. explicitly defined over all state points. The presentation will focus on examining feasible options for estimating the analysis error covariance and for defining the initial ensemble perturbations.

**Nachiketa Chakraborty**, Peter Jan van Leeuwen, Michael de Caria, and Manuel Pulido

Time varying processes in nature are often complex with non-linear and non-gaussian components. Complexity of environments and processes make it hard to disentangle different causal mechanisms which drives the observed time-series. It also makes it harder to make forecasts. The standard ways of studying causal relation in the geosciences which includes information theoretic measures of causation as well as predictive framework have deficiencies when applied to non-linear dynamical systems. Here we focus on investigating building a predictive causal framework that allows us to make predictions in simpler systems in a consistent way. We use a Bayesian framework to embed causal measures akin to mutual information from information theory to quantify relations between different random processes in this system. We examine causal relations in toy models and simple systems with a view to eventually applying to the interocean exchange problem in the Indian, the South Atlantic and the Southern Ocean.

**Maxime Conjard**and Henning Omre

The challenge in data assimilation for models representing spatio-temporal phenomena is made harder when the spatial histogram of the variable of interest appears with multiple modes. Pollution source identification constitutes one example where the pollution release represents an extreme event in a fairly homogeneous background. Consequently, our prior belief is that the spatial histogram is bimodal. The traditional Kalman model is based on a Gaussian initial distribution and Gauss-linear dynamic and observation models. This model is contained in the class of Gaussian distribution and is therefore analytically tractable. These properties that make its strenght also render it unsuitable for representing multimodality. To address the issue, we define the selection Kalman model. It is based on a selection-Gaussian initial distribution and Gauss-linear dynamic and observation models. The selection-Gaussian distribution may represent multimodality, skewness and peakedness. It can be seen as a generalization of the Gaussian distribution. The proposed selection Kalman model is contained in the class of selection-Gaussian distributions and therefore analytically tractable. The recursive algorithm used for assessing the selection Kalman model is specified. We present a synthetic case study of spatio-temporal inversion of an initial state containing an extreme event. The study is inspired by pollution monitoring. The results suggest that the use of the selection Kalman model offers significant improvements compared to the traditional Kalman model when reconstructing discontinuous initial states.

**Antoine Bernigaud**, Serge Gratton, Flavia Lenti, Ehouarn Simon, and Oumaima Sohab

We introduce a new formulation of the 4DVAR objective function by using as a penalty term a p-norm with 1 < p < 2. So far, only the 2-norm, the 1-norm or a mixed of both have been considered as regularization term. This approach is motivated by the nature of the problems encountered in data assimilation, for which such a norm may be more suited to tackle the distribution of the variables. It also aims at making a compromise between the 2-norm that tends to oversmooth the solution or produce Gibbs oscillations, and the 1-norm that tends to "oversparcify" it, in addition to making the problem non-smooth.

The performance of the proposed technique are assessed for different p-values by twin experiments on a linear advection equation. The experiments are then conducted using two different true states in order to assess the performances of the p-norm regularized 4DVAR algorithm in sparse (rectangular function) and "almost" sparse cases (rectangular function with a smoother slope). In this setup, the background and the measurements noise covariance are known.

In order to minimize the 4DVAR objective function with a p-norm as a regularization term we use a gradient descent algorithm that requires the use of duality operators to work on a non-euclidean space. Indeed, Rn together with the p-norm (1 < p < 2) is a Banach space. Finally, to tune the regularization parameter appearing in the formulation of the objective function, we use the Morozov's discrepancy principle.

**Arthur Filoche**, Julien Brajard, Anastase Charantonis, and Dominique Béréziat

The analogy between data assimilation and machine learning has already been shown and is still being investigated to address the problem of improving physics-based models. Even though both techniques learn from data, machine learning focuses on inferring model parameters while data assimilation concentrates on hidden system state estimation with the help of a dynamical model.

Also, neural networks and more precisely ResNet-like architectures can be seen as dynamical systems and numerical schemes, respectively. They are now considered state of the art in a vast amount of tasks involving spatio-temporal forecasting. But to train such networks, one needs dense and representative data which is rarely the case in earth sciences. At the same time, data assimilation offers a proper Bayesian framework allowing to learn from partial, noisy and indirect observations. Thus, each of this field can profit from the other by providing either a learnable class of dynamical models or dense data sets.

In this work, we benefit from powerful and flexible tools provided by the deep learning community based on automatic differentiation that are clearly suitable for variational data assimilation, avoiding explicit adjoint modelling. We use a hybrid model divided into 2 terms. The first term is a numerical scheme that comes from the discretisation of physics-based equations, the second is a convolutional neural network that represents the unresolved part of the dynamics. From the Data Assimilation point of view, our network can be seen as a particular parametrisation of the model error. We then jointly learn this parameterisation and estimate hidden system states within a variational data assimilation scheme. Indirectly, the issue of incorporating physical knowledge into machine learning models is also addressed.

We show that the hybrid model improves forecast skill compared to traditional data assimilation techniques. The generalisation of the method on different models and data will also be discussed.

**Haonan Ren**, Peter Jan Van Leeuwen, and Javier Amezcua

Data assimilation has been often performed under the perfect model assumption known as the strong-constraint setting. There is an increasing number of researches accounting for the model errors, the weak-constrain setting, but often with different degrees of approximation or simplification without knowing their impact on the data assimilation results. We investigate what effect inaccurate model errors, in particular, the an inaccurate time correlation, can have on data assimilation results, with a Kalman Smoother and the Ensemble Kalman Smoother.

We choose a linear auto-regressive model for the experiment. We assume the true state of the system has the correct and fixed correlation time-scale ω_{r} in the model errors, and the prior or the background generated by the model contains the model error with the fixed, guessed time-scale ω_{g} which differs from the correct one and is also used in the data assimilation process. There are 10 variables in the system and we separate the simulation period into multiple time-windows. And we use a fairly large ensemble size (up to 200 ensemble members) to improve the accuracy of the data assimilation results. In order to evaluate the performance of the EnKS with auto-correlated model errors, we calculate the ratio of root-mean-square error over the spread of all ensemble members.

The results with a single observation at the end of the simulation time-window show that, using an underestimated correlation time-scale leads to overestimated spread of the ensemble, and with an overestimated time-scale, the results show underestimation in the ensemble spread. However, with very dense observation frequency, observing every time-step for instance, the results are completely opposite to the results with a single observation. In order to understand the results, we derive the expression for the true posterior state covariance and the posterior covariance using the incorrect decorrelation time-scale. We do this for a Kalman Smoother to avoid the sampling uncertainties. The results are richer than expected and highly dependent on the observation frequency. From the analytical solution of the analysis, we find that the RMSE is a function of both ω_{r}_{ }and ω_{g}, and the spread or the variance only depends on ω_{g}. We also find that the analyzed variance is not always a monotonically increasing function of ω_{g}, and it also depends on the observation frequency. In general, the results show the effect of the correlated model error and the incorrect correlation time-scale on data assimilation result, which is also affected by the observation frequency.

**Jalisha Theanutti Kallingal**, Marko Scholze, Janne Rinne, and Johan Lindstrom

Wetlands in the boreal zone are a significant source of atmospheric methane, and hence they have been intensively studied with mechanistic models for the assessment of methane dynamics. The arctic-enabled dynamic global vegetation model LPJ-GUESS is one of the models that allow quantification and understanding of the natural methane fluxes at various scales ranging from local to regional and global, but with several uncertainties. Complexity in the underlying environmental processes, warming driven alternative paths of meteorological phenomena and changes in hydrological and vegetation conditions are exigent for a calibrated and optimised LPJ-GUESS. In this study, we used the Markov chain Monte Carlo (using Metropolis-Hastings formula) algorithm to quantify the uncertainties of LPJ-GUESS. Application of this method allows greater search of the posterior distribution, leading to a more complete characterisation of the posterior distribution with reduced risk of sample impoverishment. We will present first results from an assimilation experiment optimising LPJ-GUESS model process parameters using the flux measurement data from 2005 to 2015 from the Siikaneva wetlands in southern Finland. We analyse the parameter efficiency of LPJ-GUESS by looking into the posterior parameter distributions, parameter correlations, and the interconnections of the processes they control. As a part of this work, knowledge about how the methane data can constrain the parameters and processes is derived.

**Martin Verlaan**, Xiaohui Wang, and Hai Xiang Lin

Previous development of a parameter estimation scheme for a Global Tide and Surge Model (GTSM) showed that accurate estimation of the parameters is currently limited by the memory use of the analysis step and the computational demand. Because the estimation algorithm solver requires storage of the model output matching each observation for each parameter (or ensemble member), the requirement of memory storage gets out of control as the model simulation time increases, the model output and observation matrix become too large. The popular approach of localization does not work here because the tides propagate all over the globe in days, while parameter estimation requires weeks at least. Proper Orthogonal Decomposition (POD) is a useful technique to reduce the high dimension system with a smaller linear subspace. Singular values decomposition (SVD) is one of the methods to derive the POD modes, which is generally applied for space patterns. In this study, we focus on the application of POD in time patterns by using SVD to reduce the dimension in time patterns. As expected, the time patterns show a strong resemblance to the tidal constituents, but the same method is likely to work for a wider range of problems, which indicate that the memory requirements can be reduced dramatically by projection the model output and observations onto the time-POD patterns.

**Jeffrey Anderson**, Nancy Collins, Moha El Gharamti, Timothy Hoar, Kevin Raeder, Frederic Castruccio, Jingjing LIang, John Lin, James McCreight, Seongjin Noh, Brett Raczka, and Arezoo Arezoo Rfieeinasab

The Data Assimilation Research Testbed (DART) is a community facility for ensemble data assimilation developed and maintained by the National Center for Atmospheric Research (NCAR). DART provides ensemble data assimilation capabilities for NCAR community earth system models and many other prediction models. It is straightforward to add interfaces for new models and new observations to DART.

DART provides traditional ensemble data assimilation algorithms that implicitly assume Gaussianity and linearity. Traditional algorithms can still work when these assumptions are violated. However, it is possible to greatly improve results by extending ensemble algorithms to explicitly account for aspects of nonlinearity and non-Gaussianity. Two new algorithms have been added to DART. 1). Anamorphosis transforms variables to make the assimilation problem more linear and Gaussian before transforming posterior estimates back to the original model variables; 2). The marginal correction rank histogram filter (MCRHF) directly represents arbitrary non-Gaussian distributions. These methods are particularly valuable for data assimilation for bounded quantities like tracers or streamflow.

DART is being applied to a number of novel applications. Examples in the poster include 1). An eddy-resolving global ocean ensemble reanalysis with the POP ocean model and an ensemble optimal interpolation; 2). The WRF-Hydro/DART system now includes a multi-parametric ensemble, anamorphosis, and spatially-correlated noise for the forcing fields. 3). Results from the Carbon Monitoring System over Mountains using CLM5 to assimilate remotely-sensed observations (LAI, biomass, and SIF) for a field site in Colorado; 4). Assimilation of MODIS snow cover fraction and daily GRACE total water storage data and its impact on soil moisture using the DART/NOAH-MP system. 5). An ensemble atmospheric reanalysis using the CAM general circulation model.

**Mariusz Pagowski**, Cory Martin, Bo Huang, Daryl Kleist, and Shobha Kondragunta

In 2016 NOAA chose the FV3 (Finite Volume) dynamical core as a basis for its future global modeling system. For aerosol modeling this dynamical core was supplemented with GFS (Global Forecast System) physics and coupled through an interface with GOCART (Goddard Global Ozone Chemistry Aerosol Radiation and Transport) parameterization. The assimilation methodology relies on a hybrid variational-ensemble approach within the newly developed model-agnostic JEDI (Joint Effort for Data assimilation Integration) framework. Observations include 550 nm AOD retrievals from VIIRS (Visible Infrared Imaging Radiometer Suite) instruments on polar-orbiting SNPP and NOAA-20 satellites. The system is under development and early its results are compared with NASA'a MERRA-2 and ECMWF's CAMSiRA reanalyses.

**Youmin Tang**

In this study, the predictability of the Madden-Julian Oscillation (MJO) is investigated using the coupled Community Earth System Model (CESM) and the climatically relevant singular vector (CSV) method. The CSV method is an ensemble-based strategy to calculate the optimal growth of the initial error on the climate scale. We focus on the CSV analysis of MJO initialized at phase II, facilitating the investigation of the effect of the initial errors of the sea surface temperature (SST) in the Indian Ocean on it. Six different MJO events are chosen as the study cases to ensure the robustness of the results.

The results indicate that for all the study cases, the optimal perturbation structure of the SST, denoted by the leading mode of the singular vectors (SVs), is a meridional dipole-like pattern between the Bay of Bengal and the southern central Indian Ocean. The MJO signal tends to be more converged and significant in the Eastern Hemisphere while the model is perturbed by leading SV. The moist static energy analysis results indicate that the eastward propagation is much more evident in the terms of vertical advection and radiation flux than others. Therefore, the SV perturbation can strengthen and converge the MJO signal mostly by increasing the vertical advection of the moist static energy.

Further, the sensitivity studies indicate that the structure of the leading SV is not sensitive to the initial states, which suggests that we might not need to calculate SVs for each initial time in constructing the ensemble prediction, significantly saving computational time in the operational forecast systems.

**Ondřej Tichý**and Václav Šmídl

**y**= M

**x**+

**e**, where

**y**is the measurement vector which is typically in the form of gamma dose rates or concentrations, M is the source-receptor-sensitivity (SRS) matrix,

**x**is the unknown source term to be estimated, and

**e**is the model residue. The SRS matrix M is computed using an atmospheric transport model coupled with meteorological reanalyses. The inverse problem is typically ill-conditioned due to number of uncertainties, hence, the estimation of the source term is not straightforward and additional information, e.g. in the form of regularization or the prior source term, is often needed. Besides, traditional techniques rely on assumption that the SRS matrix is correct which is not realistic due to the number of approximations made during its computation. Therefore, we propose relaxation of the inverse model using introduction of the term Δ

_{M}such as

**y**= ( M+ Δ

_{M})

**x**+

**e**leading to non-linear inverse problem formulation, where Δ

_{M}can be, as an example, parametric perturbation of the SRS matrix M in the spatial or temporal domain. We estimate parameters of this perturbation together with solving the inverse problem using variational Bayes procedure. The method will be validated on synthetic dataset as well as demonstrated on real case scenario such as the controlled tracer experiment ETEX or episode of ruthenium-106 release over the Europe on the fall of 2017.

**Andres Yarce**, Santiago Lopez, Diego Acosta, Olga Lucia Quintero, Nicolas Pinel, Arjo Segers, and Arnold Heemink

Chemical Transport Models (CTMs) simulate the emission, transformation, and transport of atmospheric chemical species, providing concentration and deposition estimates. While greatly sophisticated, these are still imperfect representations of reality. Data Assimilation (DA), a technique whereby observations are integrated into the simulations, helps alleviate the models' weaknesses, improving their simulation outputs and enabling parameter and state estimation. The variational DA method is an efficient approach for large-scale parameter and state estimation, but it is not straightforward to implement due to the need for a tangent linear matrix of the adjoint model forecast operator. To circumvent this difficulty, the ensemble-based 4DEnVar DA technique was used in this work.

Daily NO_{2} observations from the TROPOspheric Monitoring Instrument (TROPOMI) at resolutions of 3x5 km were acquired for 2019 and assimilated into the LOTOS-EUROS CTM. Due to the scarcity of ground-based monitoring stations for atmospheric gases in Colombia, especially outside urban areas, satellite data provide an attractive alternative for DA.

The 4DEnVar DA was first evaluated via the Design of Experiments (DOE) methodology with the Lorenz96 model assimilating synthetic data. Different parameters were changed (ensemble number, spread, forcing factor and width of the assimilation time window) according to a complete 2^{4} factorial design followed by a Box Behnken design, providing an empirical model that guided the selection about how to modify those tuning parameters. The evaluation criteria used to test the 4DEnVar DA performance was the Root-Mean-Square (RMS) error between the analysis step and the synthetic data. Once this methodology was implemented, it was scaled up to the high-dimensional LOTOS-EUROS experiment.

The setup for the LOTOS-EUROS DA experiment was simplified in terms of domain area, chemical species of interest, dominant dynamics and considerations about how to perturb the parameters or initial conditions. A range of ensemble-members generated from perturbed parameters or input initial states were studied in conjunction with ensemble inflation experiments and Singular Value Decomposition projections, characterizing the degeneracy of the Gaussian assumption through the time propagation of the ensemble. Additionally, a complimentary analysis of this Gaussian ensemble degeneration was performed using the Shapiro-Wilk and Kolmogorov-Smirnov normality tests, which permitted a rational selection of the spin-up time of the model before the start of the assimilation window and the DA window size.

The assimilation of satellite NO_{2} observations into LOTOS-EUROS made possible the estimation of parameters and states. Before the DA, the non-assimilated model overestimated the magnitude of the observation, this technique improves the simulation in the sense that the analysis result approaches the observations reducing the RMS. Through this methodology, it was possible to circumvent the absence of an adjoint model associated with the chemical components of this CTM. To our knowledge, this is the first application of ensemble variational DA on a CTM for the Northwestern South America region.

**Alison Fowler**, Jozef Skákala, and Stefano Ciavatta

Monitoring biogeochemistry in shelf seas is of great significance for the economy, ecosystems understanding and climate studies. Data assimilation can aid the realism of marine biogeochemistry models by incorporating information from observations. An important source of information about phytoplankton groups and total chlorophyll is available from the ESA OC-CCI (ocean colour - climate change initiative) dataset.

For any assimilation system to be successful it is important to accurately represent all sources of data uncertainty. For the ocean colour product, the propagation of errors throughout the ocean colour algorithm makes the characterisation of the uncertainty challenging. However, the problem can be simplified by assuming that the uncertainty is a function of optical water type (OWT), which characterises the water column of each observed pixel in terms of their reflectance properties.

Within this work we apply the well-known Desroziers et al. (2005) consistency diagnostics to the Met Office’s NEMOVAR 3D-VAR DA system used to create daily biogeochemistry forecasts on the North-West European Shelf. The derived estimates of monthly ocean colour error covariances stratified by OWT are compared to previously derived estimates of the root mean square errors and biases using in-situ data match ups (Brewin et al. 2017). It is found that the agreement between the two estimates of the error variances have a strong seasonal and OWT dependence. The error correlations (which can only be estimated with the Desroziers’ method) in some instances are found to be significant out to a few 100km particularly for more turbid waters during the spring bloom. The reliability and limitation of these two estimates of the ocean colour uncertainty are discussed along with the implications for the future assimilation of ocean colour products and for ecosystem and climate studies.

**William Crawford**, Sergey Frolov, Justin McLay, Carolyn Reynolds, Craig Bishop, Benjamin Ruston, and Neil Barton

The presented work will illustrate the impact of analysis correction based additive inflation (ACAI) on atmospheric forecasts. ACAI uses analysis corrections from the NAVGEM data assimilation system as a representation of model error and is shown to simultaneously improve ensemble spread-skill, reduce model bias and improve the RMS error in the ensemble mean. Results are presented from a myriad of experiments exercising ACAI in stand-alone NAVGEM forecasts using two different ensemble systems; (1) the current operational EPS at FNMOC based on the ensemble transform method and (2) the Navy-ESPC EPS based on perturbed observations. The method of relaxation-to-prior-perturbations (RTPP) has also been implemented in the Navy-ESPC EPS and is shown to further improve the ensemble spread-skill relationship by allowing variance generated during the forecast to impact the initial-time ensemble variance in the subsequent cycle. Results from a simplified implementation of ACAI in the NAVGEM deterministic system will also be shown and indicate positive impact to model biases and RMSE.

**Yvonne Ruckstuhl**and Tijana Janjic

We investigate the feasibility of addressing model error by perturbing and estimating uncertain static model parameters using the localized ensemble transform Kalman filter. In particular we use the augmented state approach, where parameters are updated by observations via their correlation with observed state variables. This online approach offers a flexible, yet consistent way to better fit model variables affected by the chosen parameters to observations, while ensuring feasible model states. We show in a nearly-operational convection-permitting configuration that the prediction of clouds and precipitation with the COSMO-DE model is improved if the two dimensional roughness length parameter is estimated with the augmented state approach. Here, the targeted model error is the roughness length itself and the surface fluxes, which influence the initiation of convection. At analysis time, Gaussian noise with a specified correlation matrix is added to the roughness length to regulate the parameter spread. In the northern part of the COSMO-DE domain, where the terrain is mostly flat and assimilated surface wind measurements are dense, estimating the roughness length led to improved forecasts of up to six hours of clouds and precipitation. In the southern part of the domain, the parameter estimation was detrimental unless the correlation length scale of the Gaussian noise that is added to the roughness length is increased. The impact of the parameter estimation was found to be larger when synoptic forcing is weak and the model output is more sensitive to the roughness length.

**Julien Brajard**, Alberto Carrassi, and Laurent Bertino

The reconstruction from observations of the dynamics of high-dimensional chaotic models such as geophysical fluids is hampered by (i) the inevitably partial and noisy observations that can realistically be obtained, (ii) the need and difficulty to learn from long time series of data, and (iii) the unstable nature of the dynamics. To achieve such inference from the observations over long time series, it has recently been suggested to combine data assimilation and machine learning in several ways. We first rigorously show how to unify these approaches from a Bayesian perspective, yielding a non-trivial loss function.

Existing techniques to optimize the loss function (or simplified variants thereof) are re-interpreted here as coordinate descent schemes. The expectation-maximization (EM) method is used to estimate jointly the most likely model and model error statistics. The main algorithm alternates two steps: first, a posterior ensemble is derived using a traditional data assimilation step using an ensemble Kalman smoother (EnKS); second, both the surrogate model and the model error are updated using machine learning tools, a quasi-Newton optimizer, and analytical formula. In our case, the spatially extended surrogate model is formalized as a neural network with convolutional layers leveraging on the locality of the dynamics.

This scheme has been successfully tested on two low-order chaotic models with distinct identifiability, namely the 40-variable and the two-scale Lorenz models. Additionally, an approximate algorithm is tested to mitigate the numerical cost, yielding similar performances. Using indicators that probe short-term and asymptotic properties of the surrogate model, we investigate the sensitivity of the inference to the length of the training window, to the observation error magnitude, to the density of the monitoring network, and to the lag of the EnKS. In these iterative schemes, model error statistics are automatically adjusted to the improvement of the surrogate model dynamics. The outcome of the minimization is not only a deterministic surrogate model but also its associated stochastic correction, representative of the uncertainty attached to the deterministic part and which accounts for residual model errors.

**Jose M Gonzalez-Ondina**, Lewis Sampson, and Georgy Shapiro

Current operational ocean modelling systems often use variational data assimilation (DA) to improve the skill of the ocean predictions by combining the numerical model with observational data. Many modern methods are derivatives of objective (optimal) interpolation techniques developed by L. S. Gandin in the 1950s, which requires computation of the background error covariance matrix (BECM), and much research has been devoted into overcoming the difficulties surrounding its calculation and improving its accuracy. In practice, due to time and memory constraints, the BECM is never fully computed. Instead, a simplified model is used, where the correlation at each point is modelled using a simple function while the variance and length scales are computed using error estimation methods such as the Hollingsworth-Lonnberg or the NMC (National Meteorological Centre). Usually, the correlation is assumed to be horizontally isotropic, or to have a predefined anisotropy based on latitude. However, observations indicate that horizontal diffusion is sometimes anisotropic, hence this has to be propagated into BECM. It is suggested that including these anisotropies would improve the accuracy of the model predictions.

We present a new method to compute the BECM which allows to extract horizontal anisotropic components from observational data. Our method, unlike current techniques, is fundamentally multidimensional and can be applied to 2D or 3D sets of un-binned data. It also works better than other methods when observations are sparse, so there is no penalty when trying to extract the additional anisotropic components from the data.

Data Assimilation tools like NEMOVar use a matrix decomposition technique for the BECM in order to minimise the cost function. Our method is well suited to work with this type of decomposition, producing the different components of the decomposition which can be readily used by NEMOVar.

We have been able to show the spatial stability of our method to quantify anisotropy in areas of sparse observations. While also demonstrating the importance of including anisotropic representation within the background error. Using the coastal regions of the Arabian Sea, it is possible to analyse where improvements to diffusion can be included. Further extensions of this method could lead to a fully anisotropic diffusion operator for the calculation of BECM in NEMOVar. However further testing and optimization are needed to correctly implement this into operational assimilation systems.

**Ziqing Zu**, Xueming Zhu, and Hui Wang

Based on ROMS and Ensemble Optimal Interpolation (EnOI) method, the South China Sea operational Oceanography Forecasting System (SCSOFS) is implemented in National Marine Environmental Forecasting Center (NMEFC), to provide the forecast of the currents, temperature and salinity in South China Sea for the future 5 days. Recently, a systematic modification has been carried out to SCSOFS to improve its forecast skill.

For the data assimilation system, new methods have been implemented, such as using Increment Analysis Update (IAU) and First Guess at Appropriate Time (FGAT), using a high-pass filter to evaluate the background error, assimilating multi-source observations, using non-uniform localization radius. In addition, the respective contribution of each method will also be discussed.

An optimization system is implemented for evaluating the values of physical parameters in ROMS, to remove the long-term bias of simulation. Argo temperature profiles is assimilated in the first half of 2017, to obtain the optimal coefficients of horizontal/vertical viscosity/diffusion and linear bottom drag. An independent validation from July of 2017 to December of 2018 shows that the simulation is improved using the optimal values.

**Youmin Tang**and Yaling Wu

In this study, we developed a flow-dependent ensemble-based targeted observation method, by minimizing the analysis error variance under the framework of Ensemble Kalman filter (EnKF) data assimilation system. This method estimates the background error statistics as a flow dependent function. The covariance localization is also introduced for computing efficiency and alleviating the spurious observations. As a test bed, an optimal observation array of sea level anomalies (SLA) is designed for its seasonal prediction over the tropical Indian Ocean (TIO) region. Furthermore, the observing system simulation experiments (OSSEs) is used to verify the resultant optimal observational array using our recently developed coupled data assimilation system. A comparison between this flow-dependent method and the traditional method is also given.

### NP5.2 – New approaches to predictions and predictability estimation for geophysical fluid

**Bedartha Goswami**, Adam Hartland, Chaoyong Hu, Sebastian Hoepker, Bethany R. S. Fox, Norbert Marwan, and Sebastian F. M. Breitenbach

The concentration of trace elements such as Ni, Co, and Cu in a stalagmite is determined by (i) the amount of these elements present in so-called organic-metal complexes (OMCs) that trap the ionic forms of such elements in the dripwater, and (ii) the amount that is able to decay from the OMCs into the aqueous phase, from where the elements can adsorb to the growing stalagmite surface (and remain captured within the stalagmite crystal structure). A statistical treatment of the decay of a population of trace element ions from OMCs allow us to model the rates at which the dripwater dropped from the roof of the cave on to the stalagmite’s surface. The problem is however made challenging due to: (i) the lack of reliable monitoring data that quantifies the relationship between OMC trace metal ion concentration and stalagmite trace metal ion concentration, and (ii) the presence of chronological uncertainties in our estimates of trace element concentrations at past time points from the depth-based measurements along the stalagmite. We present here a semi-heuristic, semi-theoretical approach that estimates dripwater rates using a theoretical model based on the population-level chemical kinetics of trace element decay from OMCs, and a heuristic choice of calibration data sets based on precipitation and temperature from nearby weather station data. Our approach is applied to trace metal data from the Heshang Cave in southeastern China, and we are able to reconstruct a driprate proxy time series — a first quantitative hydrological proxy record presented along with well-defined estimates of uncertainty.

**Fei Zheng**, Jin-Yi Yu, and Jiang Zhu

The tropical Pacific has experienced a new type of El Niño, which has occurred particularly frequently during the last decade and is referred to as the central Pacific (CP) El Niño. Various coupled models with different degrees of complexities have been used to make real-time El Niño predictions, but large uncertainties still exist in the forecasts. It is still not yet known how much of the uncertainty is specifically related to the new CP type of El Niño and how much is common to both this type and the conventional Eastern Pacific (EP) type of El Niño. In this study, the deterministic performance of an El Niño-Southern Oscillation (ENSO) ensemble prediction system (EPS) is examined for these two types of El Niño. Ensemble hindcasts are performed for the nine EP El Niño events and twelve CP El Niño events that have occurred since 1950. The results show that (1) the skill scores for the EP events are significantly better than those for the CP events at all lead times; (2) the systematic forecast biases come mostly from the prediction of the CP events; and (3) the systematic error is characterized by an overly warm eastern Pacific during the spring season, indicating a stronger spring prediction barrier for the CP El Niño. Further improvements of coupled atmosphere-ocean models in CP El Niño prediction should be recognized as a major challenge and high-priority task for the climate prediction community.

**Tobias Braun**, Norbert Marwan, Vishnu R. Unni, Raman I. Sujith, and Juergen Kurths

We propose Lacunarity as a novel recurrence quantification measure and apply it in the context of dynamical regime transitions. Many complex real-world systems exhibit abrupt regime shifts. We carry out a recurrence plot based analysis for different paradigmatic systems and thermoacoustic combustion time series in order to demonstrate the ability of our method to detect dynamical transitions on variable temporal scales. Lacunarity is usually interpreted as a measure of ‘gappiness’ of an arbitrary spatial pattern. In application to recurrence plots, it quantifies the degree of heterogenity in the temporal recurrent patterns. Our method succeeds to distinguish states of varying dynamical complexity in presence of noise and short time series length. In contrast to traditional recurrence quantifiers, no specification of minimal line lengths is required and features beyond the scope of line structures can be accounted for. Applied to acoustic pressure fluctuation time series, it captures both the rich variability in dynamical complexity and detects shifts of characteristic time scales.

**Dmitry Mukhin**and Abdel Hannachi

We suggest a method for nonlinear analysis of atmospheric circulation regimes in the middle latitudes. The method is based on the kernel principal component analysis allowing to separate principal modes of dynamics entangled in data. We propose a new kernel function accounting specifics of large-scale wave patterns in the mid-latitude atmosphere. First, capabilities of the method are shown by the analysis of the 3-layer quasi-geostrophic model of the Northern hemisphere atmosphere: a statistically significant set of modes can be detected by the method from relatively short (several thousand days) time series. Next, we consider reanalysis data of wintertime geopotential height anomalies over the Northern hemisphere from 1950 to the present. The principal components obtained uncover several recurrent and persistent wave structures which are associated with different weather regimes. We find that there is a pronounced inter-annual and decadal variability in the dominance of different modes in different years. Possible climatic and external forcings which impact such variability as well as long-term predictability of anomalous weather seasons based on the obtained components are discussed.

**Balasubramanya Nadiga**

Whether it is turbulence fluid flows or climate variability, there is a big gap between our ability to develop understanding of underlying phenomena/processes and our ability to produce skillful predictions. We focus on near-term prediction of climate as an example. In this context, the state-of-the-art is such that we are able to predict how 30-year global averages of surface temperature will change, but we are unable to predict shorter time scale regional changes. We investigate a range of deep learning approaches to the problem ranging from reservoir computing to deep convolutional Long Short-Term Memory network architectures. The best performing architectures are seen to be capable of predicting an Earth System Model’s leading modes of global temperature variability with prediction lead times of up to a year. This approach is proposed as a useful practical tool for climate prediction. Further insight into the difficulty of the prediction problem is provided by considering the Lorenz '63 model: Long prediction horizons seen when the system is fully observed is seen to be progressively degraded as the system is less thoroughly observed, while noting the difficulty of fully observing the earth system

**Lingjiang Tao**, Wansuo Duan, and Stephane Vannitsem

Observations indicate that there exist two types of El Niño events: one is the EP-El Niño with a warming center in the eastern tropical Pacific, and the other is the CP-El Niño with large positive SST anomalies in the central tropical Pacific. Most current numerical models show low skills in identifying the El Niño diversity. The present study examines the dynamical properties of the ENSO forecast system NFSV-ICM which combines an intermediate complexity ENSO model (ICM) with a nonlinear forcing singular vector (NFSV)-tendency perturbation forecast model. This system is able to distinguish different types of El Niño in simulations and predictions. It is shown that the NFSV-ICM system is able to capture the horizontal distribution of the SST anomalies and their amplitudes in the mature phase of not only EP-El Niño but also CP-El Niño. At the same time, the NFSV-ICM is able to describe the evolution of SST anomalies associated with the two types of El Niño up to at least two-season lead time, while the corresponding forecasts with the ICM is only limited to at most one-season lead time. These improvements are associated with the modifications of atmospheric and ocean processes described by the ICM through the NFSV-tendency perturbations. In particular, the thermocline and zonal advection feedback are strongly modified and improve the conditions of emergence of both the EP- and CP-El Niño events. The NFSV-ICM therefore provides a useful platform for studying ENSO dynamics and predictability associated with El Niño diversities.

**Shraddha Gupta**, Jürgen Kurths, and Florian Pappenberger

Every point on the Earth’s surface is a dynamical system which behaves in a complex way while interacting with other dynamical systems. Network theory captures this feature of climate to study the collective behaviour of these interacting systems giving new insights into the problem. Recently, climate networks have been a promising approach to the study of climate phenomena such as El Niño, Indian monsoon, etc. These phenomena, however, occur over a long period of time. Weather phenomena such as tropical cyclones (TCs) that are relatively short-lived, destructive events are a major concern to life and property especially for densely populated coastlines such as in the North Indian Ocean (NIO) basin. Here, we study TCs in the NIO basin by constructing climate networks using the ERA5 Sea Surface Temperature and Air temperature at 1000 hPa. We analyze these networks using the percolation framework for the post-monsoon (October-November-December) season which experiences a high frequency of TCs every year. We find significant signatures of TCs in the network structure which appear as abrupt discontinuities in the percolation-based parameters during the period of a TC. This shows the potential of climate networks towards forecasting of tropical cyclones.

This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 813844.

**Kun Liu**, Jingyi Liu, Huiqin Hu, Wuhong Guo, and Baolong Cui

The sensitive area of targeted observation for the short-term prediction of the vertical thermal structure in the summer Yellow Sea is investigated by utilizing the Conditional Nonlinear Optimal Perturbation (CNOP) method and a adjoint-free algorithm with the Regional Ocean Modeling System. We use a vertical integration scheme of temperature to locate the sensitive area, in which reducing the initial errors are expected to yield great improvements in vertical thermal structure prediction of the verification area. We perform a series of sensitivity experiments to evaluate the effectiveness of the identified sensitive area. Our results show that, initially adding random perturbations in the sensitive area have the greatest negative effects on the prediction than in other areas (eg. the verification area, regions east and northeast of the verification area). Moreover, Observing System Simulation Experiments (OSSEs) indicate that, eliminating the initial errors in the sensitive area can lead to a more refined prediction than in other selected areas (including the verification area itself). Our study suggests that implementing targeted observation is a feasible way to improve the short-term prediction of the vertical thermal structure in the summer Yellow Sea.

**Stéphane Vannitsem**and Wansuo Duan

The use of coupled Backward Lyapunov vectors (BLv) for ensemble forecast is demonstrated in a coupled ocean-atmosphere system of reduced order, the Modular Arbitrary Order Ocean-Atmosphere sytem (MAOOAM). It is found that the best set of BLvs to build a coupled ocean-atmosphere forecasting system are the ones associated with near-neutral or slightly negative Lyapunov exponents. This counter intuitive result is related to the fact that these sets display larger projections on the ocean variables than the others, leading to an appropriate spread for the ocean, and at the same time a rapid transfer of these errors toward the most unstable BLvs affecting predominantly the atmosphere is experienced. The latter dynamics is a natural property of any generic perturbation in nonlinear chaotic dynamical systems, allowing for a reliable spread with the atmosphere too. The implications of these results for operational ensemble forecasts in coupled ocean-atmosphere systems are discussed.

**Wansuo Duan**, and Hui Xu

The present study uses the nonlinear singular vector (NFSV) approach to identify the optimally-growing tendency perturbations of the Weather Research and Forecasting (WRF) model for tropical cyclone (TC) intensity forecasts. For nine selected TC cases, the NFSV-tendency perturbations of the WRF model, including components of potential temperature and/or moisture, are calculated when TC intensities are forecasted with a 24-hour lead time, and their respective potential temperature components are demonstrated to have more impact on the TC intensity forecasts. The perturbations coherently show barotropic structure around the central location of the TCs at the 24-hour lead time, and their dominant energies concentrate in the middle layers of the atmosphere. Moreover, such structures do not depend on TC intensities and subsequent development of the TC. The NFSV-tendency perturbations may indicate that the model uncertainty that is represented by tendency perturbations but associated with the inner-core of TCs, makes larger contributions to the TC intensity forecast uncertainty. Further analysis shows that the TC intensity forecast skill could be greatly improved as preferentially superimposing an appropriate tendency perturbation associated with the sensitivity of NFSVs to correct the model, even if using a WRF with coarse resolution.

**Meiyi Hou**and Xiefei Zhi

Different types of El Niño-Southern Oscillation (ENSO) predictions are sensitive to the initial errors in different key areas in the Pacific Ocean. And it is known that the prediction can be improved by removing the initial errors by using assimilation methods. However yet, few studies have quantified to what extent can different types of ENSO predictions be improved by assimilating variable in different key areas. In Hou et.al (2019), 4 types of ocean temperature initial error patterns were classified for two types of El Niño prediction. It was indicated that initial errors in the north Pacific, covering the Victoria Mode region, along with south Pacific, covering the South Pacific Meridional Mode region, and subsurface layer of western equatorial Pacific have strong influence on the ENSO prediction. Following the data analysis method and the initial error patterns they proposed, we assimilate ocean temperature in these three key areas of Pacific Ocean by using CMIP5 pi-control dataset and particle filter method. Most EP- and CP-El Niño predictions in December are improved after assimilating the ocean temperature in southeast Pacific, north Pacific and western equatorial Pacific from January to March. Specially, for the prediction ensemble which contains EP(CP)-type-1 initial errors, the EP(CP)-El Niño prediction skill raises the most after assimilating the Tropical Pacific temperature, comparing with the result of assimilating the south Pacific and north Pacific. As for the prediction ensemble which contains EP-type-2 initial errors, which present similar pattern to EP-type-1 but with opposite sign, the EP-El Niño prediction skill increases the most by assimilating the north Pacific temperature. The results verify that the initial errors in the north Pacific exert contrary influences on the ENSO prediction with that in the southeast Pacific and western tropical Pacific. In addition, the initial errors in the north Pacific is more of concern for the SST prediction in the central tropical Pacific in December, while those in the southeast Pacific and tropical western Pacific is more related to the SST prediction in the central-eastern tropical Pacific. In conclusion, to better predict the types of El Niño, attentions should be paid to the initial ocean temperature accuracy not only in the tropical Pacific but also in the north and south Pacific.

**Niklas Boers**, Bedartha Goswami, Aljoscha Rheinwalt, Bodo Bookhagen, Brian Hoskins, and Jürgen Kurths

Extreme rainfall events are often coupled across long spatial distances due to atmospheric teleconnections. Revealing such linkages is important for our understanding of extreme events and related atmospheric circulation patterns, but also for enhancing the forecast skill of such events [1]. Here, we present recent results [2] on how complex networks can be employed to discover extreme rainfall teleconnections from high-resolution satellite data. Our method allows to quantitatively distinguish regional weather systems from global-scale teleconnections coupling the individual weather systems. Several lines of evidence suggest that the most relevant mechanisms for global-scale teleconnections of extreme rainfall events are related to atmospheric Rossby waves. We exemplify our approach with a focus on extreme rainfall events in the mountain regions of South-Central Asia (including Northern Pakistan and India), and show that they are statistically significantly coupled to preceding events in Europe as well as succeeding events in eastern China. An analysis of the corresponding atmospheric circulation patterns shows that a previously revealed, quasi-stationary Rossby wave termed the ‘silk road pattern’ [3] is responsible for this instance of long-range coupling between extreme rainfall events. Overall, our findings give new insights into the connections between atmospheric Rossby waves and extreme rainfall events, and thus into the potential predictability of related natural hazards. Moreover, they give promising clues in how to constrain state-of-the-art climate models with respect to their simulation of extreme rainfall.

[1] B. Hoskins: The potential for skill across the range of the seamless weather-climate prediction problem: A stimulus for our science, QJRMS 2013

[2] N. Boers, B. Goswami, A. Rheinwalt, B. Bookhagen, B. Hoskins, J. Kurths: Complex networks reveal global pattern of extreme-rainfall teleconnections, Nature 2019

[3] T. Enomoto, B. Hoskins, Y. Matsuda: The formation mechanism of the Bonin high in August, QJRMS 2003

**Liang Shi**, Ruiqiang Ding, and Yu-heng Tseng

The skills of most ENSO prediction models have declined significantly since 2000. This decline may be due to a weakening of the correlation between tropical predictors and ENSO. Moreover, the effects of extratropical ocean variability on ENSO have increased during this period. To improve ENSO predictability, we investigate the influence of the tropical-extratropical Atlantic and Pacific sea surface temperature(SST) on ENSO during the periods of pre-2000 and post-2000. We find that the influence of the northern tropical Atlantic(NTA) SST on ENSO has significantly increase since 2000. Meanwhile, there is a much earlier and stronger SST responses between NTA SST and ENSO over the central-eastern Pacific during June–July–August in the post-2000 period compared with the pre-2000 period. Furthermore, the extratropical Pacific SST predictors for ENSO still retain a ~10-month lead time after 2000. We use SST signals in the extratropical Atlantic and Pacific to predict ENSO using a statistical prediction model. These results reveal a significant improvement in ENSO prediction skills. These results indicate that the Atlantic and Pacific SSTAs can make substantial contributions to ENSO prediction, and can be further used to enhance ENSO predictability after 2000.

**Junjie Ma**and Wansuo Duan

The optimal perturbation method is a beneficial way to generate ensemble members to be used in ensemble forecasting. With orthogonal optimal perturbation, orthogonal conditional nonlinear optimal perturbations (O-CNOPs) generating initial perturbations and orthogonal nonlinear forcing singular vectors (O-NFSVs) generating model perturbations are two kinds of skillful ensemble forecasting methods. There is main disadvantage that O-CNOPs and O-NFSVs generate optimal perturbation members may need a lot of time, but in practical weather prediction, the ensemble members usually need to be generated quickly. In order to benefit from O-CNOPs and O-NFSVs, as well as considering the cost of calculation, therefore, we present a way with the big data and machine learning thinking to simplify the process of the optimal perturbation ensemble methods. Using the historical samples and their optimal perturbations to establish a database, we look for the historical sample which is analogous to what need to be forecasted currently from the database by using the convolutional neural network (CNN). In comparison with using optimization algorithm to get O-CNOPs and O-NFSVs directly, this way gets O-CNOPs and O-NFSVs faster which still obtain acceptable prediction performance. In addition, once the CNN model is trained completely, the cost of time for prediction will be saved. We illustrate the advantage by numerical simulations of a Lorenz 96 model.

Further more, based on above study, some comparison of the ensemble forecasting skill of O-CNOPs and O-NFSVs has been done, and there are three results for the reference: (1) in the early stage (1-6 days), the O-CNOPs method perform more skillfully, and in the later stage (6-12 days), the O-NFSVs method perform more skillfully; (2) within 1-5 days, if the development of analysis error is bigger than or close to the average value of the analysis error development of historical samples, the O-CNOPs method is preferred, else the O-NFSVs method is preferred; (3) within 0-3 days, if the development of energy is bigger than or close to the average value of the energy development of the historical samples, the O-CNOPs method is preferred, else the O-NFVS method is preferred. Next, further work is required to examine and explore more and deeper research using machine learning in ensemble forecasting studies of atmosphere and other systems.

**Andrey Gavrilov**, Sergey Kravtsov, Dmitry Mukhin, Evgeny Loskutov, and Alexander Feigin

According to recent study [1], the current state-of-the-art climate models lack the substantial part of internal multidecadal climate signal which is observed in the 20th century surface air temperature reanalysis data as a global stadium wave (GSW). In the presented work we further investigate this phenomenon using the recently developed method [2] of empirical spatio-temporal data decomposition into linear dynamical modes (LDMs). The important property of LDMs is their ability to take into account the time scales of the system evolution (they are extracted from observed dataset by the Bayesian optimization technique) better than some other linear techniques, e.g. traditional empirical orthogonal function decomposition. Like any linear decomposition, it provides the time series of principal components and corresponding spatial patterns.

We modify the initially developed LDM decomposition to make it possible to take into account a prescribed external forcing (like CO2 emissions, sun activity etc.) and then find part of variability which may be considered as an internal climate dynamics decomposed into set of modes with different time scales, and hence may be helpful in GSW interpretation. The results of applying the method to the 20th century surface air temperature with different ways of forcing inclusion will be presented and discussed.

1. Kravtsov, S., Grimm, C., & Gu, S. (2018). Global-scale multidecadal variability missing in state-of-the-art climate models. Npj Climate and Atmospheric Science, 1(1), 34. https://doi.org/10.1038/s41612-018-0044-6

2. Gavrilov, A., Seleznev, A., Mukhin, D., Loskutov, E., Feigin, A., & Kurths, J. (2018). Linear dynamical modes as new variables for data-driven ENSO forecast. Climate Dynamics. https://doi.org/10.1007/s00382-018-4255-7

**Qian Zhou**, Yunfei Zhang, Junya Hu, and Wansuo Duan

Considering the effects of initial uncertainty on the ENSO forecast, ensemble forecasts method is applied in the latest version of ENSO forecast system in National Marine Environmental Forecasting Center (NMEFC, China). The currently operational ENSO forecasts system of NMEFC is established based on the CESM model, with initialization and data assimilation.

First, leading five Singular Vectors (SV) are obtained using the climatological SST empirical singular vector method, and a SV based ensemble forecasts system is . However, the SVs can only present the initial errors that have the fasted error growth rates in a linear assumption, while ENSO and its forecasting system both are nonlinear. So, Conditional Nonlinear Optimal Perturbations (CNOP), which is has the largest error growth at the prediction time in a nonlinear scenario, is used to replace the leading SV, while other 4 SVs are kept to construct a CNOP-SV based ensemble forecast system. The hindcasts of ENSO from 1982 to 2017 shows that, the ENSO prediction skills of both SV based and CNOP-SV based ENSO ensemble forecasts are improved when compared with the old forecasting system, moreover, the CNOP-SV based ensemble forecast system has a much larger spread, showing higher prediction skills.

**Junya Hu**, Wansuo Duan, and Qian Zhou

The “spring predictability barrier” (SPB) is a well-known characteristic of ENSO prediction, which has been widely studied for El Niño events. However, due to the nonlinearity of the coupled ocean–atmosphere system and the asymmetries between El Niño and La Niña, it is worthy to investigate the SPB for La Niña events and reveal their differences with El Niño. This study investigates the season-dependent predictability of sea surface temperature (SST) for La Niña events by exploring initial error growth in a perfect model scenario within the Community Earth System Model. The results show that for the prediction through the spring season, the prediction errors caused by initial errors have a season-dependent evolution and induce an SPB for La Niña events. Two types of initial errors that often yield the SPB phenomenon are identified: the first are type-1 initial errors showing positive SST errors in the central-eastern equatorial Pacific accompanied by a large positive error in the upper layers of the eastern equatorial Pacific. The second are type-2 errors presenting an SST pattern with positive errors in the southeastern equatorial Pacific and a west–east dipole pattern in the subsurface ocean. The type-1 errors exhibit an evolving mode similar to the growth phase of an El Niño-like event, while the type-2 initially experience a La Niña-like decay and then a transition to the growth phase of an El Niño-like event. Both types of initial errors cause positive prediction errors for Niño3 SST and under-predict the corresponding La Niña events. The resultant prediction errors of type-1 errors are owing to the growth of the initial errors in the upper layers of the eastern equatorial Pacific. For the type-2 errors, the prediction errors originate from the initial errors in the subsurface layers of the western equatorial Pacific. These two regions may represent the sensitive areas of targeted observation for La Niña prediction. In addition, the type-2 errors in the equatorial regions are enlarged by the recharge process from 10°N in the central Pacific during the eastward propagation. Therefore, the off-equatorial regions around 10°N in the central Pacific may represent another sensitive area of La Niña prediction. Additional observations may be prioritized in these identified sensitive areas to better predict La Niña events.

**Evgeny Loskutov**, Valery Vdovin, Andrey Gavrilov, Dmitry Mukhin, and Alexander Feigin

We investigate the Middle Pleistocene Transition (MPT) - a rapid change in the periodicity of the Pleistocene glacial cycles from 41 kyr to about 100 kyr, which occurred about a million years ago - using the data-driven model [1]. Here we estimate stability of the model using a novel concept of interval stability [2-4], referring to the behavior of the perturbed model during a finite time interval. In a few words we define the class of 'safe' perturbations after which the system (our data-driven model) returns back to the initial dynamical regime and 'unsafe' perturbation of minimal amplitude needed to disrupt the system.

We demonstrate that the MPT is likely associated with decreasing of the climate system's interval stability to rapid disturbances (millennial and shorter). This confirms the statement made in the paper [1] that the main factor in the onset of the long-period glacial cycles is strongly nonlinear oscillations induced by the short-scale variability.

- D. Mukhin, A. Gavrilov, E. Loskutov, J. Kurths, A. Feigin. Bayesian Data Analysis for Revealing Causes of the Middle Pleistocene Transition. Scientific Reports, 9 7328 (2019).
- P. Menck, J. Heitzig, N. Marwan, J. Kurths. How basin stability complements the linear-stability paradigm. Nature Phys, 9 89–92 (2013).
- V. Klinshov, V. Nekorkin, J. Kurths. Stability threshold approach for complex dynamical systems. New Journal Physics, 18 013004 (2016).
- V. Klinshov, S. Kirillov, J. Kurths, V. Nekorkin. Interval stability for complex systems. New Journal Physics, 18 013004 (2018).

**Lin Jiang**and Wansuo Duan

Previous studies show that the kinetic energy of mesoscale eddies (MEs) accounts for more than 80% of the global ocean energy. The theoretical study and numerical simulation of MEs will enable us to better understand the dynamics of ocean circulation. Weiss and Grooms (2017) found that assimilating uniform observations taken over MEs is much better than assimilating a subset of observations on a regular grid for improving prediction skill of SSH associated with ocean state. In the present study, we use a conditional nonlinear optimal perturbation (CNOP) approach to investigate the sensitivity of the ocean state sea surface height (SSH) predictions on MEs with a two-layer quasi-geostrophic model and show the optimal assimilating scheme. In the study, the CNOPs of SSH predictions are first computed. It is found that, if one regards the regions covered by the grid points with large values of CNOPs as sensitive area of SSH predictions, the sensitive areas are mainly located on MEs. Furthermore, the stronger the MEs, the more the MEs grid points covered by the sensitive area. Especially, these grid points associated with sensitive areas are not uniformly distributed over the MEs. It is obvious that the predictions of SSH are quite sensitive to the initialization of MEs (especially that of the particular region of large values of CNOPs for strong MEs, rather than of the uniformly distributed grid points over MEs). Therefore, an appropriate initialization of MEs is much helpful for improving the prediction accuracy of SSH. And the CNOPs of SSH prediction here may provide useful information on how to improve initialization of MEs.

**Panayiotis A Varotsos**

By analyzing the seismicity in the new time domain termed natural time [1], the entropy changes of seismicity before major earthquakes have been studied. It was found [2-5] that the key quantity is the entropy change ΔS under time reversal, which is minimized a few months before major earthquakes such as the M9.0 Tohoku earthquake [2] on 11 March 2011 and the M8.2 Chiapas earthquake [3] in Mexico on 7 September 2017; accompanied by an abrupt increase of its fluctuations [4,5]. Here we discuss how these fluctuations may lead to a procedure through which the occurrence time of an impending mainshock can be estimated [6].

References

1. Varotsos P.A., Sarlis N.V. and Skordas E.S., *Natural Time Analysis: The new view of time. Precursory Seismic Electric Signals, Earthquakes and other Complex Time-Series* (Springer-Verlag, Berlin Heidelberg) 2011.

2. N. V. Sarlis, E. S. Skordas, and P. A. Varotsos, "A remarkable change of the entropy of seismicity in natural time under time reversal before the super-giant M9 Tohoku earthquake on 11 March 2011", EPL (Europhysics Letters), 124 (2018), 29001.

3. N. V. Sarlis, E. S. Skordas P. A. Varotsos, A. Ramírez-Rojas, E. L. Flores-Márquez, "Natural time analysis: On the deadly Mexico M8.2 earthquake on 7 September 2017", Physica A 506 (2018), 625-634.

4. P. A. Varotsos, N. V. Sarlis and E. S. Skordas, "Tsallis Entropy Index q and the Complexity Measure of Seismicity in Natural Time under Time Reversal before the M9 Tohoku Earthquake in 2011", Entropy 20 (2018), 757.

5. A. Ramírez-Rojas, E. L. Flores-Márquez, N. V. Sarlis and P. A. Varotsos, "The Complexity Measures Associated with the Fluctuations of the Entropy in Natural Time before the Deadly México M8.2 Earthquake on 7 September 2017", Entropy 20 (2018), 477.

6. E. S. Skordas, N. V. Sarlis and P. A. Varotsos “Identifying the occurrence time of an impending major earthquake by means of the fluctuations of the entropy change under time reversal”, EPL (Europhysics Letters), *in press*.

**Alexander Feigin**, Aleksei Seleznev, Dmitry Mukhin, Andrey Gavrilov, and Evgeny Loskutov

We suggest a new method for construction of data-driven dynamical models from observed multidimensional time series. The method is based on a recurrent neural network (RNN) with specific structure, which allows for the joint reconstruction of both a low-dimensional embedding for dynamical components in the data and an operator describing the low-dimensional evolution of the system. The key link of the method is a Bayesian optimization of both model structure and the hypothesis about the data generating law, which is needed for constructing the cost function for model learning. The form of the model we propose allows us to construct a stochastic dynamical system of moderate dimension that copies dynamical properties of the original high-dimensional system. An advantage of the proposed method is the data-adaptive properties of the RNN model: it is based on the adjustable nonlinear elements and has easily scalable structure. The combination of the RNN with the Bayesian optimization procedure efficiently provides the model with statistically significant nonlinearity and dimension.

The method developed for the model optimization aims to detect the long-term connections between system’s states – the memory of the system: the cost-function used for model learning is constructed taking into account this factor. In particular, in the case of absence of interaction between the dynamical component and noise, the method provides unbiased reconstruction of the hidden deterministic system. In the opposite case when the noise has strong impact on the dynamics, the method yield a model in the form of a nonlinear stochastic map determining the Markovian process with memory. Bayesian approach used for selecting both the optimal model’s structure and the appropriate cost function allows to obtain the statistically significant inferences about the dynamical signal in data as well as its interaction with the noise components.

Data driven model derived from the relatively short time series of the QG3 model – the high dimensional nonlinear system producing chaotic behavior – is shown be able to serve as a good simulator for the QG3 LFV components. The statistically significant recurrent states of the QG3 model, i.e. the well-known teleconnections in NH, are all reproduced by the model obtained. Moreover, statistics of the residence times of the model near these states is very close to the corresponding statistics of the original QG3 model. These results demonstrate that the method can be useful in modeling the variability of the real atmosphere.

The work was supported by the Russian Science Foundation (Grant No. 19-42-04121).

**Jeong Ock Lim**, Kyung-On Boo, HaeJin Kong, Sanghee Jun, Jeong-Hyun Park, and Hyun-Suk Kang

A typhoon is a high-impact weather phenomenon that causes serious damage to people and property when landing on the Korean peninsula. Therefore, accurate prediction of typhoon intensity and track is very important in establishing damage prevention measures.

KMA has applied its own typhoon initialization process (KMA bogusing) for Global Data Assimilation and Prediction System (GDPS) to produce realistic initial fields since 2010, when the Unified Model (UM) of UK Met Office (UKMO) was introduced as an operational model. If the typhoon intensity of the background is weaker than the observation, the KMA bogussing process generates horizontally spread mean sea level pressures by using TC warning center’s advisory and use it for data assimilation.

In June 2018, GDPS has been upgraded based on PS40 N1280 of UKMO with significantly increased its horizontal resolution from 17 km to 10 km. The new version of GDPS showed improved performance in TC intensity predictions. Since the model simulates tropical cyclone intensity strong enough, we investigated the impact of typhoon initialization on the predictability of the new version GDPS.

**Jingyi Liu**, Wuhong Guo, Baolong Cui, Kun Liu, and Huiqin Hu

Targeted observation is an appealing procedure to improve oceanic model predictions by taking additional assimilation of collected measurements. However, studies on targeted observation in the oceanic field have been largely based on modeling efforts, and there is a need for field validating observations. Here, we report the preparatory work of a field campaign, which is designed based on the identified sensitive area by the Conditional Nonlinear Optimal Perturbation (CNOP) approach, to improve the short-range summer thermal structures prediction in the Yellow Sea (YS). We firstly simulated the hindcasting (2016-2018) temperature structures in the summertime, and found that the locations of the sensitive areas are generally consistent in space for each hindcast year. Then, we introduced the technique of multiple-assimilation and the definition of time-varying sensitive area, and designed observing strategies for the YS summer campaign. Observing System Simulation Experiments (OSSEs) were conducted prior to address the plan on field campaign in the Yellow Sea in August 2019. Results show that, reducing the initial errors in the sensitive area can lead to more improvement on thermal structures prediction than that in other area.

**Ville Leinonen**, Petri Tiitta, Olli Sippula, Hendryk Czech, Ari Leskinen, Juha Karvanen, Sini Isokääntä, and Santtu Mikkonen

Aerosols and their transformation process in atmosphere have significant effects on climate. Transformation process is a complex combination of physical and chemical reactions. Multiple oxidizing agents and other factors, such as radiation, affect the transformation process. Characterization of these factors and their strength is a problem, where advanced methods might help to gain more understanding.

In this work, we modeled transformation of wood combustion emission measured in the environmental chamber by using causal modeling (Pearl, 2009). The aim of the study was to use state-of-the-art causal discovery methods to search causal pathways between measured variables: precursors and particle products. The data used in the modelling are introduced in Tiitta et al. (2016).

In addition to wood combustion experiments, we simulated artificial datasets to understand abilities of the model. We wanted to evaluate the accuracy of our model to confirm the correct structure between variables and reproduce the measured transformation. This helps us to understand the model performance in real datasets.

We found that model could reproduce the measured evolution well. The structure between emission parts was not completely matching to prior assumption. Usually incorrect predictors in the modeled structure are highly correlated with correct causes.

References:

Pearl, J.: Causality : Models, Reasoning and Inference., Cambridge University Press., 2009.

Tiitta et al., Atmos. Chem. Phys., 16, 13251-13269, 2016.

**Wuhong Guo**, Baolong Cui, and Jingyi Liu

Focusing on the transfer and evolvement of initial perturbation from temperature of ocean to the underwater acoustic propagation,comparing with the remote sensing data,optimizing Ocean-Acoustic Coupled Model, the reliability of this model is verified. On this basis,global and regional error development experiments are carried out by adding perturbation on the initial temperature of a controlled test. The results show that after 5 days evolution of the initial temperature field,the global perturbation of the propagation loss is saturated, and the perturbation structure is basically consistent with the law of the dynamic ocean.For the target area in Kuroshio region, the initial perturbation in the upstream region is the fastest. This conclusion can provide the basis for the adaptive observation of ocean acoustics.

**Baolong Cui**and Wuhong Guo

Focusing on the rapid prediction of acoustic field uncertainty in environment with temporal and spatial sound speed perturbation, evolvement of sound speed structure over time is predicted based on the ocean-acoustic coupled model to obtain the uncertainty distribution of the vertical structure of sound speed. Further, a method combining the arbitrary polynomial chaos expansion with the empirical orthogonal function is proposed to reduce the dimensionality of uncertain parameters and to obtain the uncertainty distribution of the acoustic field. Simulations have shown that the computational complexity can be reduced by 2 orders of magnitude compared to the conventional polynomial chaos expansion while ensures the same precision. Moreover, the computational complexity is not influenced by the complexity of the sound speed profile. The acoustic field and uncertainty predicted in uncertain environment by proposed method also have been tested with the experimental data.

### NP5.4 – Advances in statistical post-processing, blending and verification of deterministic and ensemble forecasts

**Moritz N. Lang**, Sebastian Lerch, Georg J. Mayr, Thorsten Simon, Reto Stauffer, and Achim Zeileis

Non-homogeneous regression is a frequently-used post-processing method for increasing the predictive skill of probabilistic ensemble weather forecasts. To adjust for seasonally varying error characteristics between ensemble forecasts and corresponding observations, different time-adaptive training schemes, including the classical sliding training window, have been developed for non-homogeneous regression. This study compares three such training approaches with the sliding-window approach for the application of post-processing near-surface air temperature forecasts across Central Europe. The predictive performance is evaluated conditional on three different groups of stations located in plains, in mountain foreland, and within mountainous terrain, as well as on a specific change in the ensemble forecast system of the European Centre for Medium-Range Weather Forecasts (ECMWF) used as input for the post-processing.

The results show that time-adaptive training schemes using data over multiple years stabilize the temporal evolution of the coefficient estimates, yielding an increased predictive performance for all station types tested compared to the classical sliding-window approach based on the most recent days only. While this may not be surprising under fully stable model conditions, it is shown that "remembering the past" from multiple years of training data is typically also superior to the classical sliding-window when the ensemble prediction system is affected by certain model changes. Thus, reducing the variance of the non-homogeneous regression estimates due to increased training data appears to be more important than reducing its bias by adapting rapidly to the most current training data only.

**Matthieu Lafaysse**, Guillaume Evin, Matthieu Vernay, Joseph Bellier, Bruno Joly, Maxime Taillardat, and Michaël Zamo

Forecasting the height of new snow (HS) is essential for avalanche hazard survey, road and ski resorts management, tourism attractiveness, etc. Meteo-France operates the PEARP-S2M probabilistic forecasting system including 35 members of the PEARP Numerical Weather Prediction system, the SAFRAN downscaling tool refining the elevation resolution in mountains, and the Crocus snowpack model representing the main physical processes in the snowpack (compaction, melting, etc.). It provides better HS forecasts than direct NWP diagnostics but exhibits significant biases and underdispersion. Therefore, a post-processing is required to be able to provide automatic forecasting products of HS from this system.

For that purpose, we compare the skill of two statistical methods (Nonhomogeneous Regression with a Censored Shifted Gamma distribution and Quantile Regression Forest), two predictor datasets for training (22-year reforecast with some discrepancies with the operational system or 3-year real time forecasts similar to the operational system) and two spatial scales of post-processing (local scale or 1000 km² regional scale).

The improvement relative to the raw forecasts is similar at both spatial scales. Thus, the regional validity of post-processing does not restrict the application at points with observations. The impact of the training dataset depends on lead time and on the evaluation criteria. The long-term reforecast improves the reliability of severe snowfall but leads to overdispersion due to a discrepancy with the initial perturbations used in the operational system. Finally, thanks to a larger number of predictors, the Quantile Regression Forest allows an improvement of forecasts for specific cases when the the rain-snow transition elevation is overestimated by the raw forecasts.

These conclusions help to choose an optimal post-processing configuration for automatic forecasts of the height of new snow and encourage the atmospheric modelling teams to develop long reforecasts as homogenous as possible with the operational systems.

**Julian Steinheuer**and Petra Friederichs

**Jonas Bhend**, Christoph Spirig, Max Hürlimann, Lionel Moret, and Mark Liniger

Weather forecasts have been steadily improving in quality over the last decades. These ongoing improvements are due to advances in numerical weather prediction (NWP) and the advent of ever more powerful supercomputers that allow simulating future weather and its uncertainty with increasing resolution and using ensemble approaches. Such physics-based computer models, however, are not free of systematic errors. Statistical postprocessing can be used to calibrate NWP forecasts to further improve forecast quality and better exploit the available information. Here we present results from several explorative deep learning studies using artificial neural networks (ANN) to calibrate high resolution forecasts of temperature, precipitation, wind, and cloud cover in Switzerland. These first attempts at ANN-based postprocessing help us to understand the strengths and weaknesses of machine learning and are the basis to build more complex and comprehensive statistical models accounting for local effects in complex terrain such as the Swiss Alps. In all cases, ANN leads to significant improvements over the direct NWP output. While the improvement is comparable in magnitude with improvements achieved with conventional postprocessing approaches, ANN-based postprocessing is easier to generalize in space for a calibration of forecasts also at unobserved sites. In addition to the results of the postprocessing, we will also discuss the lessons learned so far in using machine learning for this particular problem.

**Maurice Schmeits**, Simon Veldkamp, and Kirien Whan

Current statistical post-processing methods for providing a probabilistic forecast are not capable of using full spatial patterns from the numerical weather prediction (NWP) model output. Recent developments in deep learning (notably convolutional neural networks) have made it possible to use large gridded input data sets. This could potentially be useful in statistical post-processing, since it allows us to use more spatial information.

In this study we consider wind speed forecasts for 48 hours ahead, as provided by KNMI's Harmonie-Arome model. Convolutional neural networks, fully connected neural networks and quantile regression forests are used to obtain probabilistic wind speed forecasts. Comparing these methods shows that convolutional neural networks are more skillful than the other methods, especially for medium to higher wind speeds.

**Florian Dupuy**, Olivier Mestre, and Léo Pfitzner

Cloud cover is a crucial information for many applications such as planning land observation missions from space. However, cloud cover remains a challenging variable to forecast, and Numerical Weather Prediction (NWP) models suffer from significant biases, hence justifying the use of statistical post-processing techniques. In our application, the ground truth is a gridded cloud cover product derived from satellite observations over Europe, and predictors are spatial fields of various variables produced by ARPEGE (Météo-France global NWP) at the corresponding lead time.

In this study, ARPEGE cloud cover is post-processed using a convolutional neural network (CNN). CNN is the most popular machine learning tool to deal with images. In our case, CNN allows to integrate spatial information contained in NWP outputs. We show that a simple U-Net architecture produces significant improvements over Europe. Compared to the raw ARPEGE forecasts, MAE drops from 25.1 % to 17.8 % and RMSE decreases from 37.0 % to 31.6 %. Considering specific needs for earth observation, special interest was put on forecasts with low cloud cover conditions (< 10 %). For this particular nebulosity class, we show that hit rate jumps from 40.6 to 70.7 (which is the order of magnitude of what can be achieved using classical machine learning algorithms such as random forests) while false alarm decreases from 38.2 to 29.9. This is an excellent result, since improving hit rates by means of random forests usually also results in a slight increase of false alarms.

**Sam Allen**, Chris Ferro, and Frank Kwasniok

Raw output from deterministic numerical weather prediction models is typically subject to systematic biases. Although ensemble forecasts provide invaluable information regarding the uncertainty in a prediction, they themselves often misrepresent the weather that occurs. Given their widespread use, the need for high-quality wind speed forecasts is well-documented. Several statistical approaches have therefore been proposed to recalibrate ensembles of wind speed forecasts, including a heteroscedastic censored regression approach. An extension to this method that utilises the prevailing atmospheric flow is implemented here in a quasigeostrophic simulation study and on reforecast data. It is hoped that this regime-dependent framework can alleviate errors owing to changes in the synoptic-scale atmospheric state. When the wind speed strongly depends on the underlying weather regime, the resulting forecasts have the potential to provide substantial improvements in skill upon conventional post-processing techniques. This is particularly pertinent at longer lead times, where there is more improvement to be gained upon current methods, and in weather regimes associated with wind speeds that differ greatly from climatology. In order to realise this potential, however, an accurate prediction of the future atmospheric regime is required.

**Jon Saenz**, Sheila Carreno-Madinabeitia, Ganix Esnaola, Santos J. González-Rojí, Gabriel Ibarra-Berastegi, and Alain Ulazia

A new diagram is proposed for the verification of vector quantities generated by individual or multiple models against a set of observations. It has been designed with the idea of extending the Taylor diagram to two-dimensional vector such as currents, wind velocity, or horizontal fluxes of water vapour, salinity, energy and other geophysical variables. The diagram is based on a principal component analysis of the two-dimensional structure of the mean squared error matrix between model and observations. This matrix is separated in two parts corresponding to the bias and the relative rotation of the empirical orthogonal functions of the data. We test the performance of this new diagram identifying the differences amongst a reference dataset and different model outputs using examples wind velocities, current, vertically integrated moisture transport and wave energy flux time series. An alternative setup is also proposed with an application to the time-averaged spatial field of surface wind velocity in the Northern and Southern Hemispheres according to different reanalyses and realizations of an ensemble of CMIP5 models. The examples of the use of the Sailor diagram show that it is a tool which helps identifying errors due to the bias or the orientation of the simulated vector time series or fields. An implementation of the algorithm in form of an R package (sailoR) is already publicly available from the CRAN repository, and besides the ability to plot the individual components of the error matrix, functions in the package also allow to easily retrieve the individual components of the mean squared error.

**Swati Singh**, Kaustubh Salvi, Subimal Ghosh, and Subhankar Karmakar

The downscaling approaches: Statistical and Dynamic, developed for regional climate predictions, have both advantages and limitations. The statistical downscaling is computationally inexpensive but suffers from the violation of the assumption of stationarity in statistical (predictor-predictand) relationship. The dynamical downscaling is assumed to take care of stationarity but suffers from the biases associated with various sources. Here we propose a joint approach of both the methods by applying statistical methods: bias correction & statistical downscaling to **Coordinated Regional Climate Downscaling Experiment (**CORDEX) evaluation runs. The evaluation runs are considered as perfect simulations of CORDEX Regional Climate Models (RCMs) with the boundary conditions by ERA-Interim reanalysis data. The statistical methods are also applied to ERA-Interim reanalysis data and compared with observation data for Indian Summer Monsoon characteristics. We evaluate the ability of statistical methods under the non-stationary environment by taking the difference of years close to extreme future runs (RCP8.5) as warmer years and preindustrial runs as cooler years. We find statistical downscaling of CORDEX evaluation runs shows skill in reproducing the signal of non-stationarity. The study can be extended methods by applying statistical downscaling to CORDEX RCMs with the CMIP5 boundary conditions.

**Zhaolu Hou**, Bin Zuo, Shaoqing Zhang, Fei Huang, Ruiqiang Ding, Wansuo Duan, and Jianping Li

Numerical forecasts always have associated errors. Analogue correction methods combine numerical simulations with statistical analyses to reduce model forecast errors. However, identifying appropriate analogues remains a challenging task. Here, we use the Local Dynamical Analog (LDA) method to locate analogues and correct model forecast errors. As an example, an ENSO model forecast error correction experiment confirms that the LDA method locates more dynamical analogues of states of interest and better corrects forecast errors than do other methods. This is because the LDA method ensures similarity of the initial states and the evolution of both states. In addition, the LDA method can be applied using a scalar time series, which reduces the complexity of the dynamical system. Model forecast error correction using the LDA method provides a new approach to correcting state-dependent model errors and can be readily integrated with other advanced models.

**Stéphane Vannitsem**

For most statistical post-processing schemes used to correct weather forecasts, changes to the forecast model induce a considerable reforcasting effort. We present a new approach based on response theory to cope with slight model change. In this framework, the model change is seen as a perturbation of the original forecast model. The response theory allows then to evaluate the variation induced on the averages involved in the statistical post-processing, provided that the magnitude of this perturbation is not too large.

This approach is studied in the context of a simple quasi-geostrophic model. It provides a proof-of-concept of the potential performances of response theory in a chaotic system. The parameters of the statistical post-processing used - an *Error-in-Variables* Model Output Statistics (EVMOS) - are appropriately corrected when facing a model change. The potential application in a more operational environment is also discussed.

**Juwon Kim**, Hae-Jin Kong, and Hyuncheol Shin

Multi-model ensemble using statistical post-processing is one of the methods to provide the impact of uncertainties of the Numerical Weather Prediction (NWP) models, with low cost and better accuracy for extreme weather forecasts. Extreme weather events such as heat/cold waves, windstorms, and heavy rainfall result in severe damage in human life and properties. However, the performance of the NWP models, particularly, heavy rain forecast is still low due to the intermittent and non-Gaussian properties. The light rain tends to be overestimated and the strong rain tends to be underestimated averagely on the NWP models. Thus the multi-model ensemble using statistical post-processing is activated to correct the discrepancies between the observation and the model intensity of precipitation.

The aim of this study is to provide the improvement of precipitation forecasts in probabilistic and deterministic aspects using a multi-model ensemble method with more weights on the less error and without any bias correction. Six types of models, namely, Local Data assimilation and Prediction System (LDPS), Local ENsemble System (LENS), Global Data assimilation and Prediction System (GDPS), Ensemble Prediction System-Global (EPSG) of Korea Meteorological Administration (KMA), the single and ensemble models of European Centre for Medium-Range Weather Forecasts (ECMWF), are used to blend. The preliminary results of the multi-model ensemble show similar results to the ECMWF ensemble mean in deterministic for 3-hourly accumulated precipitation over the East Asia and the middle of the performance among individual models in probabilistic over the South Korea. More details of the methodology, results, and improvements will be discussed in the presentation.

**Iris Odak Plenkovic**, Suzana Panezic, and Endi Keresturi

Even the state-of-the-art mesoscale models exhibit noteworthy errors, especially in the complex terrain. Therefore, it is useful to include post-processing methods in the forecasting system to further reduce starting model errors at locations where measurements are available.

The analog-based method (ABM) is a point-based post-processing approach which consists of two steps. The first step is to find the most similar past numerical weather predictions (analogs) over several variables (predictors) and the second is to form an analog ensemble (AnEn) out of the corresponding observations.

The ABM is thoroughly tested using the wind speed NWP input, focusing on the complex terrain. Since August 2019 it is used in a test operational mode at Croatian Meteorological and Hydrological Service. The setup includes 15 members wind speed, wind gusts and temperature ensemble predictions for approximately 50 stations using the 2-year training dataset. The preliminary results show that the ABM implementation is successful, reducing the error and improving the skill of the raw model. Additionally, it is found that the ABM predictions of wind speed and gusts optimally need more predictors than the temperature predictions. Finally, the forecasting system shows the best result in the coastal region for the temperature predictions, while the best results for the wind speed are achieved in the nearly flat continental terrain situated more inland.

**Alexander Pasternack**, and Henning Rust

In weather and climate science statistical modeling is applied for manifold problems. Due to the increasing number of input variables, overfitting can easily deteriorate the performance for model predictions. In order to avoid this, it is often meaningful to apply model selection approaches. Since conventional approaches can be very time-consuming especially for many predictors, we are using the boosting approach, which combines model selection and parameter estimation. This iterative algorithm identifies and updates in each step only the most important coefficient, such that in the end most important predictor variables have non-zero coefficients and less relevant variables are ignored.

Boosting has been originally developed for classification problems but has also been extended and used for other applications; i.a. non-homogeneous gaussian regression. Based on the non-homogeneous boosting proposed by Messner et al. (2016), which is used to model mean and variance of a forecast distribution simultaneously, we have developed a boosting algorithm for a non-stationary Generalized Extreme Value distribution (GEV). Thus, it is possible to identify the most relevant predictor variables for location, scale and shape parameter concurrently. We apply this algorithm to various toy model simulations to assess the effect of this novel approach.

**Maxime Taillardat**and Olivier Mestre

Statistical post-processing of ensemble forecasts, from simple linear regressions to more sophisticated techniques, is now a well-known procedure in order to correct biased and misdispersed ensemble weather predictions. However, practical applications in National Weather Services is still in its infancy compared to deterministic post-processing. This paper presents two different applications of ensemble post-processing using machine learning at an industrial scale. The first is a station-based post-processing of surface temperature in a medium resolution ensemble system. The second is a gridded post-processing of hourly rainfall amounts in a high resolution ensemble prediction system. The techniques used rely on quantile regression forests (QRF) and ensemble copula coupling (ECC), chosen for their robustness and simplicity of training whatever the variable subject to calibration.

Moreover, some variants of classical techniques used such as QRF or ECC have been developed in order to adjust to operational constraints. A forecast anomaly-based QRF is used for temperature for a better prediction of cold and heat waves. A variant of ECC for hourly rainfall is built, accounting for more realistic longer rainfall accumulations. It is shown that forecast quality as well as forecast value is improved compared to the raw ensemble. At last, comments about model size and computation time are made.