Detecting regime transitions of the nocturnal and Polar near-surface temperature inversion
Amandine Kaiser, Davide Faranda, Sebastian Krumscheid, Danijel Belušić, Nikki Vercauteren
DDetecting regime transitions of the nocturnal and Polarnear-surface temperature inversion
Amandine Kaiser , Davide Faranda , Sebastian Krumscheid , Danijel Belu˘si´c , and NikkiVercauteren FB Mathematik und Informatik, Freie Universit¨at Berlin, 14195 Berlin, Germany LSCE-IPSL, CEA Saclay l’Orme des Merisiers, CNRS UMR 8212 CEA-CNRS-UVSQ,Universit´e Paris-Saclay, 91191 Gif-sur-Yvette, France London Mathematical Laboratory, London, United Kingdom Department of Mathematics, RWTH Aachen, Aachen, Germany SMHI, Norrk¨oping, SwedenMarch 24, 2020
Abstract
Many natural systems undergo critical transitions, i.e. sudden shifts from one dynamical regime toanother. In the climate system, the atmospheric boundary layer can experience sudden transitions be-tween fully turbulent states and quiescent, quasi-laminar states. Such rapid transitions are observed inPolar regions or at night when the atmospheric boundary layer is stably stratified, and they have im-portant consequences in the strength of mixing with the higher levels of the atmosphere. To analyze thestable boundary layer, many approaches rely on the identification of regimes that are commonly denotedas weakly and very stable regimes. Detecting transitions between the regimes is crucial for modelingpurposes.In this work a combination of methods from dynamical systems and statistical modeling is appliedto study these regime transitions and to develop an early-warning signal that can be applied to non-stationary field data. The presented metric aims at detecting nearing transitions by statistically quan-tifying the deviation from the dynamics expected when the system is close to a stable equilibrium. Anidealized stochastic model of near-surface inversions is used to evaluate the potential of the metric as anindicator of regime transitions. In this stochastic system, small-scale perturbations can be amplified dueto the nonlinearity, resulting in transitions between two possible equilibria of the temperature inversion.The simulations show such noise-induced regime transitions, successfully identified by the indicator. Theindicator is further applied to time series data from nocturnal and Polar meteorological measurements.
The atmospheric boundary layer (ABL) is the lowest part of the atmosphere that is directly influencedby the Earth’s surface and across which turbulent exchanges of momentum, heat and matter between thesurface and the free atmosphere occur. During daytime, surface warming leads to an unstable or convectiveboundary layer. During clear-sky nights, radiative cooling leads to a surface that is cooler than the airaloft and the ABL becomes stably stratified. The stable stratification can also arise when warm air isadvected over a colder surface, which is a frequent event in Polar regions. Turbulence in the resulting stableboundary layer (SBL) is subject to buoyant damping and is only maintained through mechanical productionof turbulent kinetic energy (TKE). Understanding and modeling the SBL is essential for regional and globalatmospheric models, yet there are many well-documented challenges to simulate stably stratified atmosphericflows [Sandu et al., 2013, Holtslag et al., 2013, LeMone et al., 2018]. One of the challenges is to developan accurate understanding and representation of distinct regimes of the SBL and transitions between them[Baas et al., 2017]. 1 a r X i v : . [ phy s i c s . a o - ph ] M a r umerous observational and modeling studies show that the SBL can be classified, to a first approxima-tion, in a weakly stable regime in which turbulence is continuous, and a very stable regime with patchy andintermittent turbulence, requiring a different modeling approach [Mahrt, 2014]. The weakly stable regimetypically occurs when cloud cover limits nocturnal radiative cooling at the land surface, or with strong windsassociated to wind shear that produces enough TKE to sustain turbulence. The vertical mixing is thereforemaintained and a well-defined boundary layer usually exists in which turbulent quantities decrease upwardsfrom the surface layer following the classical model of Monin-Obukhov similarity theory and related existingconcepts [Mahrt, 2014]. The associated temperature stratification, or temperature inversion, is weak. Thestrongly stable regime occurs with strong stratification and weak winds and does not follow the traditionalconcept of a boundary layer. Transitions from weakly stable to strongly stable regimes are caused by astrong net radiative cooling at the surface which increases the inversion strength and eventually leads tosuppressed vertical exchanges, unless winds are strong enough to maintain turbulence [van de Wiel et al.,2007]. The reduced vertical mixing results in a decoupling from the surface, such that similarity theorybreaks down [Acevedo et al., 2015]. Intermittent bursts tend to be responsible for most of the turbulenttransport [Acevedo et al., 2006, Vercauteren et al., 2016]. Such bursts alter the temperature inversion andcan sometimes drive transitions from strongly to weakly stable boundary layers.Transitions between the different SBL regimes have been found by modeling studies to be dynamicallyunpredictable. Based on a numerical model representing the exchanges between the surface and the SBLusing a very simple two-layer scheme, McNider [1995] showed the existence of bi-stable equilibria of thesystem, which can thus transition between very different states under the influence of random perturbations.Interacting nonlinear processes that lead to this bi-stability partly involve thermal processes at the landsurface, as was highlighted following the hypothesis that continuous turbulence requires the turbulence heatflux to balance the surface energy demand resulting from radiative cooling [van de Wiel et al., 2007, 2012a,2017]. According to this maximum sustainable heat flux hypothesis, a radiative heat loss that is strongerthan the maximum turbulent heat flux that can be supported by the flow with a given wind profile will leadto the cessation of turbulence [van de Wiel et al., 2012b] and thus to a regime transition. This concept is usedby van Hooijdonk et al. [2015] to show that the shear over a layer of certain thickness can predict SBL regimeswhen sufficient averaging of data is considered. Based on observations, Sun et al. [2012] identify a heightand site-dependent wind speed threshold that triggers a transition between a regime in which turbulenceincreases slowly with increasing wind speed from a regime where turbulence increases rapidly with the windspeed. The change of the relationship between the turbulence and the mean wind speed occurs abruptlyat the transition. Sun et al. [2016] attribute this difference to the turbulent energy partitioning betweenturbulent kinetic energy and turbulent potential energy (TPE): in the very stable regime, shear-inducedturbulence will have to enhance the TPE in order to counter the stable stratification before enhancing theTKE.The combined importance of the wind speed and of the surface thermal processes has also been evidencedby numerical studies using idealized single-column models of the atmosphere. Single-column models with afirst-order turbulence closure scheme [Baas et al., 2017, 2019, Holdsworth and Monahan, 2019] or a second-order closure scheme [Maroneze et al., 2019] are able to representatively simulate transitions from weaklyto strongly-stable regimes. Yet, direct numerical simulations show that transitions from strongly to weakly-stable regimes can occur following a localized, random perturbation of the flow [Donda et al., 2015]. Fieldstudies have also highlighted examples of transitions induced by small-scale perturbations of the flow [Sunet al., 2012]. In fact, a statistical classification scheme introduced by Vercauteren and Klein [2015] showsthat the SBL flow transitions between periods of strong and weak influence of small-scale, non-turbulent flowmotions on TKE production in the SBL. Such submeso-scale fluctuations of the flow (e.g. induced by variouskind of surface heterogeneity) are typically not represented in models but are important in strongly stableregimes [Vercauteren et al., 2019], and may trigger regime transitions. Stochastic modeling approaches area promising framework to analyze their impact on regime transitions. A related statistical classification ofthe Reynolds-averaged boundary-layer states introduced by Monahan et al. [2015] highlighted that regimetransitions are a common feature of SBL dynamics around the globe [Abraham and Monahan, 2019]. Regimetransitions typically take an abrupt character [Baas et al., 2019, Acevedo et al., 2019]. Predicting thetransition point remains a challenge [van Hooijdonk et al., 2016].Abrupt or critical transitions are ubiquitous in complex natural and social systems. The concept ofcritical transition is formally defined in dynamical systems theory and relates to the notion of bifurcation2Kuznetsov, 2013]. When the dynamics is controlled by a system of equations depending on an externalparameter (often called forcing), the stability of the equilibrium solutions can change abruptly and this isalso reflected on macroscopic observables of the system. Sometimes, one can have early-warning signals of atransition because the systems experience some influences of the bifurcated state before actually reaching it.The motion of a particle undergoing random fluctuations in an asymmetric double-well energy potential V is a minimal system to detect early warnings, in which each well or local minimum of the energy potentialcorresponds to a stable equilibrium state of the system. For a small fixed level of noise, the control parameteris ∆ V , the depth of the well in which the particle is located, leading to an energy barrier that the particlehas to overcome in order to transition to the second stable equilibrium. If ∆ V is large, the distributionof the positions of the particle will be quasi Gaussian and the autocorrelation function of the position ofthe particle will have an exponential decay. Conversely, if ∆ V is reduced when the particle approaches thebifurcation point, then the particle position’s distribution starts to feel the effect of the other state and thedistribution will be skewed towards the new state. Similarly, the excursions from the equilibrium positionwill become larger, increasing the autocorrelation time. Early-warning signals can then, e.g., be definedbased on the changes of the autocorrelation time.These first early-warning signs have been successfully applied to several systems with excellent re-sults [Scheffer et al., 2009, 2012], including the present context of SBL regime transitions [van Hooijdonket al., 2016]. However, sometimes, transitions can happen without detectable early-warnings [Hastings andWysham, 2010]. The main limitation of early-warning signals based on the increase of autocorrelation isthat their activation does not always correspond to a bifurcation. Indeed, if a single well potential widens, asit can occur in non-stationary systems, the distribution of a particle’s position experiences the same increasein skewness and autocorrelation function without the need of approaching a bifurcation [Lenton et al., 2012,Faranda et al., 2014]. In our context of SBL flows, non-stationarity of the energy potential governing thedynamics can be due to changes in the mean wind speed or cloud cover, for example. For these reasons,[Faranda et al., 2014] have introduced a new class of early-warning indicators based on defining a distancefrom the dynamics expected from a particle evolving in a single-well potential. The suggested indicatorstatistically quantifies the dynamical stability of the observables and was already used by Nevo et al. [2017]to show that strongly stable flow regimes are dynamically unstable and may require high-order turbulenceclosure schemes to represent the dynamics. Alternative new early warnings are based on the combination ofstatistical properties of observables when approaching the bifurcation [Chen et al., 2012].In the present analysis, we investigate if the early-warning indicator introduced by Faranda et al. [2014]can be used to detect nearing transitions between SBL flow regimes, based on both simulated data andfield measurements. We show that the conceptual model that was recently suggested by van de Wiel et al.[2017] to understand SBL regime transitions in terms of thermal coupling of the land surface is equivalentto a dynamical system representing the evolution of the temperature inversion evolving in a double-wellenergy potential. We extend this conceptual model to a stochastic model where added noise representsthe effect of natural fluctuations of the temperature inversion’s rate of change. The resilience of equilibriaof the non-random model to perturbations as well as the bifurcation points are known analytically (as wasdiscussed in van de Wiel et al. [2017]), and we thus use the simulated data to test our indicator. Additionally,the indicator relies on calculating statistical properties of the data with a moving window approach and issensitive to the choice of the window length. We suggest two complementary, data-driven but physicallyjustified approaches to define an appropriate window length for which results can be trusted. Finally, theindicator is applied to nocturnal temperature inversion data from a site in Dumosa, Australia as well as fromtemperature inversion data from Dome C, Antartica. The goal of our study is to investigate if a statistical early-warning indicator of regime transitions canbe successfully used to detect nearing regime transitions in the SBL. In section 2.1, the conceptual modelintroduced by van de Wiel et al. [2017] to study regime transitions will be introduced, along with its dynamicalstability properties. In section 2.2, the model is extended to a stochastic model in which noise representsfluctuations in the dynamics of the near-surface temperature inversion. In section 2.3, we present a statisticalindicator that was introduced in Faranda et al. [2014] and applied to SBL turbulence data in Nevo et al.32017] to estimate the dynamical equilibrium properties of time series, based on a combination of dynamicalsystems concepts and stochastic processes tools. The conceptual model describes the evolution of the near-surface temperature inversion and is used to produce time series of controlled data for which the theoreticalequilibrium properties are known.
A conceptual model was introduced by van de Wiel et al. [2017] to study regime transitions of near-surfacetemperature inversions in the nocturnal and polar atmospheric boundary layer. The authors were able todetermine a connection between the dynamical stability of the temperature inversion and the ambient windspeed U through their model and measurements. Mathematically speaking, the model is a dynamical sys-tem represented by a first order ordinary differential equation, abbreviated ODE, which describes the timeevolution of the difference between the temperature at a reference height T r and the surface temperature T s .Although the equilibrium properties of the system and the dynamical stability properties (i.e. the resilienceto perturbations) of all equilibria states were thoroughly discussed in van de Wiel et al. [2017], for the sakeof completeness we briefly introduce the model and summarize the linear stability analysis of equilibriumpoints of the resulting ODE for different values of a bifurcation parameter. The bifurcation parameter isrelated to the ambient wind speed.Assuming that the wind speed and temperature are constant at a given height z r , the following equa-tion describes the evolution of the near-surface inversion strength, based on a simple energy balance at theground surface: c v d ∆ Tdt = Q n − G − H. (1)In this energy balance model, c v is the heat capacity of the soil, ∆ T = T r − T s is the inversion strengthbetween the temperature at height z r and at the surface z s , Q n is the net long wave radiative flux (an energyloss at the surface that will be set as a constant), G is the soil heat flux (an energy storage term that willbe parameterized as a linear term), and H is the turbulent sensible heat flux (a non-linear energy transportterm that will be parameterized in the following).After parameterizing the fluxes, the model has the form: c v d ∆ Tdt = Q i − λ ∆ T − ρc p c D U ∆ T f ( R b ) , (2)in which Q i is the isothermal net radiation, λ is a lumped parameter representing all feedbacks from soil heatconduction and radiative cooling as a net linear effect, ρ is the density of air at constant pressure, c p is theheat capacity of air at constant pressure, c D = ( κln ( z r /z ) ) is the neutral drag coefficient with κ ≈ . z the roughness length and z r the reference height, U is the wind speed at height z r , R b = z r ( gT r ) ∆ TU is the bulk Richardson number, and f ( R b ) is the stability function used in Monin-Obukhovsimilarity theory.The lumped parameter λ corresponds to a linear term in the model as the soil is assumed to respondlinearly to the temperature inversion. Moreover, ∆ T · f ( R b ) is a non-linear term due to the non-lineardependence of turbulent diffusion on the vertical temperature gradient.Following van de Wiel et al. [2017], instead of analyzing the dynamical stability of the energy-balance model(2) itself, we will present the linear stability analysis of a simplified system that has a similar mathematicalstructure but is mathematically convenient to analyze. Using a cutoff, linear form for the stability function,i.e. f ( x ) = 1 − x and f ( x ) = 0 for x >
1, the simplified model is dx ( t ) dt = g (cid:0) x ( t ) (cid:1) , where g ( x ) = (cid:40) Q i − λx − Cx (1 − x ) for x ≤ Q i − λx for x > x ( t ) = x . Here, up to dimensional constants, x represents ∆ T . The parameter C will be treated asa bifurcation parameter for this simplified system. Similar types of stability functions are typically used in4umerical weather prediction tools, and the cutoff form facilitates the mathematical analysis of the model.Note that to be consistent with the original model, the stability function should include a dependence onboth the temperature and the wind speed via R b . Removing this dependence as it is done here changes someof the nonlinearity, however it makes the mathematical analysis very simple and the qualitative behavior ofthe system is similar to the original system - see van de Wiel et al. [2017], their Figure 8 and 10. In thatsense, the model loses some physical significance for mathematical convenience, but the qualitative nonlinearfeedback processes are maintained. For a deeper discussion of the model, its simplifications and the modelparameters, the reader is referred to its thorough presentation by van de Wiel et al. [2017].For fixed and physically meaningful values of Q i and λ , equation (3) can have either one, two, or threepossible equilibrium solutions depending on the fixed values (see illustration in van de Wiel et al. [2017],Figure 10-12, and related discussion for more details). The equilibrium solutions will be functions of theparameter C , which we will consider as a bifurcation parameter in the following. Physically, the case ofstrong thermal coupling between the surface and the atmosphere, corresponding to a large value of λ , resultsin one unique equilibrium solution whose value depends on C . In van de Wiel et al. [2017], it is hypothesizedthat such a case is representative of a grass site such as Cabauw, the Netherlands. The solution is linearlystable to perturbations, i.e. linear stability analysis shows that perturbed solutions are attracted back tothe equilibrium. The case of no coupling ( λ =0) leads to two equilibrium solutions, one of which is linearlystable and the other unstable to perturbations (i.e. perturbed solutions are repelled by the equilibrium). Aweak coupling strength, with an intermediate value of λ that could be representative of a snow surface oranother thermally insulated ground surface, results in three possible equilibrium solutions. The two extremesolutions are stable to perturbations, while the middle equilibrium solution is unstable. Perturbed solutionsaround the middle equilibrium will thus be attracted either by the upper or the lower equilibrium. Plottingthose three equilibrium solutions as a function of the bifurcation parameter C results in a back-folded curvewhich is qualitatively similar to observations of the temperature inversion shown as a function of wind speedat Dome C, Antartica; see Vignon et al. [2017]. The bifurcation diagram is shown in Figure 1 for parametervalues such that λ > Q i > λ , resulting in the case with three possible equilibrium solutions. By C x e stableunstable Figure 1: Bifurcation plot for simplified model. The lowest, middle and upper branches correspond respec-tively to the equilibria x e , x e and x e convention, the unstable equilibrium branch is denoted by a dashed line. In the following, we will analyzetransitions between the two stable equilibria. If the system undergoes random perturbations in this bi-stablecontext, a perturbation could drive the system sufficiently far from a stable equilibrium state so that itcomes near the unstable equilibrium and finally gets attracted by the second equilibrium. The three possibleequilibria are denoted as x e , x e and x e . In order to study such regime transitions induced by randomperturbations, the conceptual model is extended with a noise term in the following section.5 .2 Extending the conceptual model by randomization The conceptual model (3) can be equivalently written in terms of a gradient system, in which the temperatureinversion represented here by x evolves according to the influence of an underlying potential V ( x ). Therandomized model to be introduced will be based on this gradient structure. Specifically, the initial valueproblem (3) can be written as dxdt = − dVdx , x ( t ) = x , where it is easy to see that the potential is given by V ( x ) = (cid:40) x ( λ + C ) − C x − Q i x for x ≤ , λx − Q i x + C for x > . (4)The linear stability analysis discussed in the previous section can thus be understood in the sense that thetemperature inversion x equilibrates at a local minimum of a potential V . That is, an equilibrium point x e satisfies V (cid:48) ( x e ) = 0. Figure 2 sketches the form of the potential with the exemplary parameter values λ = 2 , Q i = 2 . , C = 6 .
4. Note that V ( x ) is a double-well potential in that case where each local minimumcorresponds to one of the stable equilibrium points x e and x e , while the local maximum corresponds tothe system’s unstable equilibrium x e .While the conceptual model (3) has proven very insightful to explain observed sharp transitions intemperature inversions, it only allows for regime transitions when drastic changes in the model parameters (i.e., bifurcations) occur. That is, the model is overly idealized and in reality one can expect regime transitionsto also take place due to small natural fluctuations of the temperature inversion itself in certain cases, e.g.,when the potential barrier separating the two local minima and corresponding stable equilibria is shallow.Therefore we will consider an appropriate randomized variant of the model. Specifically, we consider thestochastic differential equation (SDE) model dx = − dV ( x ) dx dt + σ dB , x ( t ) = x , (5)to account for small random perturbations to the temperature inversion’s rate of change. Here, B denotesa standard Brownian motion (i.e., a stochastic process) and σ > V is as in (4). As the randomized dynamics is characterized by the same potential, alsothe equilibrium points of the non-random model (3) will describe the dominant effects of the randomizedmodel’s dynamics. However, due to the presence of the noise, the stable equilibria of the non-randommodel (3) are not limiting points for the stochastic counter-part in (5), in the sense that the temperatureinversion may still fluctuate after reaching a stable equilibrium. The reason is that in a context of twostable equilibria (i.e. for parameter values such that the model (3) exhibits two stable equilibria, denotedearlier as the thermally weakly-coupled state), the random perturbations can trigger transitions from onestable equilibrium to another one. We will therefore refer to the formerly stable states as: metastable.Note that depending on the coupling strength and noise intensity, the likelihood of regime transitions canchange drastically and the system may or may not exhibit metastable states. The type of noise (additiveor multiplicative for example, or noise with a Levy distribution) will also affect regime transitions. In oursubsequent simulations and analyses, we will focus on the case of two metastable states with additive noiseand leave other cases for future research.The effect of these random perturbations to a metastable equilibrium point x e can be understood througha localized approximation of the original dynamics. More precisely, consider a second-order Taylor approxi-mation of the potential around an equilibrium point x e , yielding the quadratic approximate potential ˜ V : V ( x ) ≈ ˜ V ( x ) := V ( x e ) + 12 d Vdx ( x ) (cid:12)(cid:12)(cid:12)(cid:12) x = x e ( x − x e ) . For the same parameter values that were used to plot the original potential in Figure 2, the red line in thesame figure shows the approximate quadratic potential around the equilibrium value x e . Using the locallyquadratic potential, we can thus define a locally approximate dynamics for the temperature inversion byreplacing V in (5) by ˜ V , resulting in dX = − k ( X − x e ) dt + σ dB , X ( t ) = x e , .0 0.5 1.0 1.5 2.0 - . . . . x V ( x ) stableunstable ˜ V ( x ) Figure 2: Example of a potential V ( x ) (eq. 4) and its local approximation through the quadratic potential˜ V ( x ) (see text for details).where k := d Vdx ( x ) | x = x e ∈ R and X is introduced to describe the approximate dynamics of the former x .This approximate dynamics is an example of the well-studied Ornstein–Uhlenbeck process and it provides anaccurate description of the full dynamics in the neighborhood to the equilibrium point x e . Discretizing theOrnstein–Uhlenbeck process X using the Euler–Maruyama scheme with a step-size ∆ t := TL for some positiveinteger L we furthermore find that the process at discrete times t ∈ { , . . . , L } approximately satisfies thedifference equation X t = X t − − k ( X t − − x e )∆ t + σ (cid:16) B ( t ∆ t ) − B (cid:0) ( t − t (cid:1)(cid:17) , X = e e , in the sense that X t ≈ X ( t ∆ t ). By defining µ := kx e ∈ R , φ := (1 − k ∆ t ) and w t := σ (cid:2) B ( t ∆ t ) − B (cid:0) ( t − t (cid:1)(cid:3) this can be written as X t = µ + φX t − + w t , which is a so-called autoregressive model of order one, denoted AR(1), thanks to the properties of the(scaled) Brownian increments w t . Consequently, we see that the discretized Ornstein–Uhlenbeck process canbe accurately approximated by an AR(1) process. This derivation can also be found in Thomson [1987].Combining this with the observation that the Ornstein–Uhlenbeck process offered an accurate approximationto the original dynamics in the vicinity of a stable equilibrium, we can thus conclude that the local dynamicsin the neighborhood of a metastable state can be approximately described by an AR(1) process. In section 2.1 we discussed the simplified model by van de Wiel et al. [2017] which was developed to un-derstand regime transitions in near-surface temperature inversions. This model provides a hypothesis thatexplains the existence of two possible equilibria of the temperature inversion for a given wind speed. Inagreement with the randomized conceptual model introduced in the previous section, we say that a systemexhibiting at least two metastable equilibria is called metastable. In this section the goal is to describe amethodology for statistically detecting critical transitions based on time series data. For the detection weapply an indicator for the dynamical stability (i.e. the resilience to perturbations) of time series, which wasdefined by Faranda et al. [2014] and applied to SBL turbulence data in Nevo et al. [2017]. The indicator usesa combination of methods from dynamical systems and from statistical modeling. In its definition, deviationsfrom AR(1) processes in the space of so-called autoregressive-moving-average (ARMA) models are used to7uantify the dynamical stability of a time series. A time series x t , t ∈ Z , is an ARMA(p,q) process if it isstationary and can be written as x t = ν + p (cid:88) i =1 φ i x t − i + w t + q (cid:88) j =1 θ j w t − j , (6)with constant ν , coefficients φ p , θ q and { w t } being white noise with positive variance σ . The coefficients φ p and θ q additionally have to satisfy some constraints, see Brockwell and Davis [2016]. Notice that AR(1)= ARMA(1,0). Intuitively the parameters p and q are related to the memory lag of the process. The longerthe system takes to return to the equilibrium after a perturbation, the more memory we expect to observein the process. Examples of simple systems along with their ARMA(p,q) characteristics can be found inFaranda et al. [2014].In section 2.2, it was shown that the dynamics when the system is close to a stable equilibrium can beapproximated by a AR(1) process. We will assume that far from the transition from one dynamical regimeto another, the time series of a generic physical observable can be described by an ARMA(p,q) model witha reasonably low number of p, q parameters and coefficients. Indeed far from a transition, the system willtend to remain around an equilibrium despite random perturbations, and excursions from the equilibrium areshort. The idea behind the modeling assumption is that ARMA processes are an important parametric familyof stationary time series [Brockwell and Davis, 2016]. Their importance is due to their flexibility and theircapacity to describe many features of stationary time series. Thereby, choosing ARMA(p,q) processes formodeling the dynamics away from a stable state is a reasonable Ansatz. Close to a transition, the resilienceof the system to perturbations is weak and longer excursions from the equilibrium occur. The statisticalproperties (such as the shape and/or the persistence of the autocorrelation function) of the system changedrastically, leading to an expected increase of the value p + q [Faranda et al., 2014]. Based on this, we useARMA(p,q) models in the following to analyze the stability of a dynamical system. The dynamical stabilityindicator which will be defined next will be used to obtain indicators for detecting the system’s proximityto a regime transition.In order to quantify the local stability of a time series, we first slice the time series x t with a moving timewindow of fixed length τ . In other words, we obtain subsequences { x , . . . , x τ } , { x , . . . , x τ +1 } , . . . , { x t − τ +1 , . . . , x t } of the original time series that overlap. By slicing the original time series we obtain a sequence of shorter timeseries for which it is reasonable to suppose that they are amenable to ARMA modeling. In detail, we assumethat the subsequences are realizations of linear processes with Gaussian white noise which then implies thatthe process is stationary. We then fit an ARMA(p,q) model for every possible value of ( p, q ), with p ≤ p max and q ≤ q max , to these subsequences, where p max and q max are predefined thresholds. Afterwards we choosethe best fitting ARMA(p,q) model by choosing the one with the minimal Bayesian information criterion,BIC, [Schwarz, 1978] which is a commonly used and well-studied tool in statistical model selection. Assum-ing that we have the maximum likelihood estimator ˆ β := (ˆ ν, ˆ φ , . . . , ˆ φ p , ˆ θ , . . . , ˆ θ q ) of the fitted ARMA(p,q)model (which can be obtained using a so-called innovation algorithm, as it is, for example, implemented inthe ”forecast” R package [Hyndman et al., 2019] which is used for the analyses), the BIC is formally definedas BIC( p, q ) = − L ( ˆ β ) + ln( τ )( p + q + 1) , (7)where L ( ˆ β ) denotes the associated likelihood function evaluated at the maximum likelihood estimator ˆ β .The second term introduces a penalty for high-order models (i.e., those that contain more parameters) toavoid over-fitting.We reiterate that when the system is close to a metastable state, its dynamics can be well approximatedby an AR(1) process. When the system is approaching an unstable point separating multiple basins ofattraction, the approximation no longer holds as the underlying potential cannot be approximated by aquadratic potential anymore. The change in the shape of the potential introduces new correlations in thetime series, resulting in higher-order ARMA terms when fitting such a model to data.The definition of the stability indicator is based on this observation, in the sense that it assumes thatthe dynamics near a metastable state can be modeled by an ARMA(1,0) or AR(1)-process. Specifically, the8tability indicator is defined asΥ( p, q ; τ ) = 1 − exp (cid:18) −| BIC( p, q ) − BIC(1 , | τ (cid:19) . (8)For a stable state, the most likely statistical model is an AR(1) process and one expects that Υ = 0. Theindicator Υ gives a normalized distance between the stable state (Υ = 0) and the state in which the systemis. The limit Υ → p, q and τ inthe following discussion.The reliability of Υ highly depends on the choice of τ , the window length (which we will consider in numberof discrete observations in the following), and it relates to the characteristic timescales of the dynamics.Intuitively, the window length, when converted to its equivalent physical duration (i.e. the number of discreteobservations multiplied by the discrete sampling time), has to be shorter than the residence timescale in onebasin of attraction (i.e. the time spent in the neighborhood of an equilibrium before transitioning to anotherone) in order to satisfy the local stationarity, but large enough so that statistical model estimation is reliable.In winter at Dome C where the Polar winter results in a near absence of daily cycle, no preferred timescaleof residence around an equilibrium of the temperature inversion was observed (an equilibrium can remainfor several days), however the transition between two equilibria was observed to take place over a timescaleof the order of 10 hours [Baas et al., 2019]. For nocturnal flows, the residence timescale is tightly connectedto the daily cycle and could be of a few hours during the night, or the entire night. The transition betweentwo equilibria typically takes place over a duration of about a half hour. For reliable statistical estimation,multiple tests showed that a minimum window length of 20 discrete points is needed. With a sampling timeof 1 minute, that means that a moving window of approximately 20 or 30 minutes may be appropriate.In addition, the sampling frequency has to be fine enough to sample typical fluctuations of the dynamics.In the following analyses, we find a sampling frequency of one minute to be appropriate for that purpose.The characteristic timescale here is given by the timescale at which the system recovers from perturbations(which is estimated by linear stability analysis in the case where the model is known, see e.g. [van de Wielet al., 2017]), and the time interval between two observations should be close to or smaller than this quantityso that (small-scale) local fluctuations can be resolved. Since the characteristic timescales of the systemcannot be known analytically in many situations, for example when analyzing time series from atmosphericmodels or from field data, we suggest two data-driven approaches to select a window length: • In the first approach, the mean residence time around each metastable state as well as the meantransition time between the two states will be estimated based on a data clustering approach . Theobservations will be clustered in the metastable regimes and an intermediate, transition regime. Fromthe clustered data, the mean residence time in each cluster will be evaluated. This approach willprovide an upper bound to select the window length. • The second approach is motivated by the fact that the indicator Υ is obtained through a statisticalinference procedure through the definition of the BIC which involves fitting suitable ARMA processesto data. Specifically, a maximum likelihood approach is used, which assumes that subsequences aresampled from a normal distribution. To assess the validity of this statistical approach , a normality testwill be implemented as a criterion to select a window length for which the normality assumption isjustified and ARMA model estimation is reliable.Both approaches are applicable when the data-generating model is unknown. This is important in caseswhere data showing signs of metastability are available, but an underlying model is unknown. A summaryof the full algorithmic procedure used to calculate the statistical indicator is given in Appendix 4.3.
In the first approach, we suggest to use the K-means algorithm (Hartigan and Wong [1979], see pseudo codein Appendix 4.1) to select a window length for the analysis. In the context of analyzing transitions in the9emperature inversion, the idea is to cluster the data into three different clusters: data around each stablefixed point and data near the unstable fixed point (in other words, data covering the transition periodsbetween two metastable states). By that, the goal is to estimate the average time needed by the systemto transition between two metastable states. The mean residence time within each cluster is calculatedfrom the time series of cluster affiliation. We choose τ (recall that we consider it in number of discreteobservations and not in physical time) such that it is smaller than the minimal mean time spent in onecluster, which should ensure that subsequences remain mostly around one equilibrium. This value is denotedby τ Kmeans . For the simulated data in the following, each simulated time series will be assigned a windowlength τ Kmeans by this procedure. For the nocturnal dataset, we cluster the entire dataset once and obtaina length τ Kmeans of 22 points, corresponding to a duration of 22 minutes. For the Polar dataset, only onecontinuous time series during a Polar winter will be considered and assigned one value of τ Kmeans , namely10 points, corresponding to a duration of 100 minutes. This window length is insufficient to obtain reliablestatistical estimations.Note that this clustering approach to determining a residence timescale around an equilibrium is a crudeapproximation and suffers from many caveats: a high density of data close to a given value of the temperatureinversion may not necessarily relate to the existence of metastable equilibria, but could occur due to non-stationary dynamics or complex nonlinear effects, for example. Nevertheless, we use it as a first approachand future research may result in more reliable approaches. In the following analyses, the K-means procedurecan be interpreted as providing an upper bound for selecting a window length for the analysis and thus,combined with the following criterion, will offer an applicability criterion for our method.
The K-means clustering approach described above estimates the system’s physical timescales, but the sta-tistical properties of the process should also be considered for reliable calculations. To fit ARMA modelsreliably and to calculate the Bayesian information criterion for ARMA model selection, we need the under-lying process to follow a normal distribution. Consequently, we suppose that the subsequences are sampledfrom a normal distribution, at least for some window length τ . We then choose τ as the biggest windowlength such that this normality assumption holds (more precisely, such that the normality hypothesis cannotbe rejected). This value is denoted by τ AD . Specifically, the statistical test results in a p-value for eachsubsequence, and we choose the window length such that the median of the p-values of all subsequences isabove a threshold for which the null-hypothesis cannot be rejected. The normality test applied here is theAnderson-Darling Test [Anderson and Darling, 1952], abbreviated AD test, as it is, for example, more stablethan the Kolmogorov–Smirnov test [Stephens, 1974]. Further details details of the AD test are summarizedin Appendix 4.2. Similarly to the clustering approach, the window length τ AD for which normality cannotbe rejected is selected for each analyzed time series, and this value of τ AD is then used for all subsequences ofthe time series. Each of the simulated dataset, the nocturnal dataset and the Polar dataset will be assigneda single value of τ AD . The value of τ AD for the nocturnal dataset is 19 discrete points, hence 19 minuteswith a sampling frequency of one minute. For the Polar dataset, the value of τ AD is 43 points correspondingto a duration of 430 minutes. In this section we quantify the reliability of the stability indicator introduced in section 2.2.3. We start bytesting it on a controlled dataset generated by the simplified model for near-surface temperature inversion(see section 2.2.1) and then proceed by applying Υ, the stability indicator, to observational data. In the testswe use the auto.arima() function from the ”forecast” R package [Hyndman et al., 2019]. The auto.arima()function fits ARMA(p,q) models by calculating the maximum likelihood estimators for a given model order(using the innovation algorithm mentioned earlier). It calculates the corresponding BIC (using the definition(7)) for all ARMA(p,q) models with p ≤ p max and q ≤ q max , where p max and q max are thresholds set to10 in our application, and then it chooses the ARMA model with the minimal BIC value. This procedureis repeated for each subsequence of data, using the moving window approach, and the minimal BIC valueleads to the optimal ARMA(p,q) model to represent the given subsequence.10 .1 Simulated time series To generate the simulated data, we use the conceptual randomized model (5), which we recall here for thereader’s convenience: dx ( t ) = − dV (cid:0) x ( t ) (cid:1) dx dt + σdB ( t ) , x ( t ) = x , ≤ t ≤ T , where V ( x ) is the energy potential defined in (4). That is, the data-generating model reads dx = (cid:40) ( Q i − λx − Cx (1 − x )) dt + σdB, x ≤ , ( Q i − λx ) dt + σdB, x > . (9)The SDE model (9) is approximated path-wise (i.e., for each realization of the driving Brownian path) usingthe Euler-Maruyama scheme.For the purpose of testing the accuracy of the Υ indicator and its potential to detect nearing regimetransitions, one realization { x t } of the stochastic process is used for each fixed value of the bifurcationparameters C . Multiple fixed values of C are used, resulting in one timeseries per value of C . The initialparameters are set to t = 0 and x ( t ) = min { x e i | i = 1 , , } where x e i are the three equilibria of the system.To generate the controlled data set the model parameters are set to λ = 2, Q i = 2 . σ = 0 .
35. Thevalue of C is varied between C = 5 . C = 7 . . C . The simulations are ran for n = 2000 time steps of size ∆ t = 0 .
01. The amplitude ofthe noise, or diffusion coefficient, σ = 0 .
35 is chosen as it resulted in trajectories for which regime transitioncould be observed on the time interval [0 , T = n ∆ t = 20]. The range for C is chosen because for these valuesthe time series shows frequent transitions from one metastable state to another. To choose the windowlength τ we apply both the K-means Algorithm (section 2.2.3.2.3.1) and the Anderson Darling Test (section2.2.3.2.3.2). The length τ is determined individually for each simulation, i.e. for each fixed value of C . TheK-means algorithm can be used to estimate the amount of discrete observation points covering the transitiontime. We set the cluster number to three as we expect three equilibria. The results of the clusteringalgorithm are exemplary shown in figure 3 for C = 6 .
4. Note that t ∈ [0 , n · ∆ t ] = [0 , . . . . . t ˜ x Figure 3: Clustered simulated time series for C = 6 . τ KMeans = 94timesteps, with a discrete timestep ∆ t = 0 .
01. Hence the window duration corresponds to 0 .
94 time units.express our window length τ in number of discrete points in the following. In this case the equilibria are x e = 0 .
46 (metastable) , x e = 0 .
97 (unstable), and x e = 1 .
25 (metastable). The cluster centers, estimated11y the K-means algorithm, are 0.46, 0.97, and 1.31 which are a close approximation of the equilibria.Therefore, we expect a good estimation for the amount of points covering the transition. The average timespend in each cluster are (for C = 6 . mean ( T ) = 112 . mean ( T ) = 94, and mean ( T ) = 286 .
67, where mean ( T i ) is the average time spent without observed transitions in cluster i ∈ { , , } expressed in numberof discrete points. The minimal mean residence time is thus mean ( T ) = 94 and provides an upper boundto select a window length that respects the timescales of the system. The window length τ is thus chosensuch that it is smaller than the minimal average time spent in one cluster, i.e. for C = 6 . τ <τ KMeans := min { mean ( T i ) | i = 1 , , } = 94. For all tested C , we choose τ = min { mean ( T i ) | i = 1 , , }− τ ). By applying Υ to the data generated by the simplified model with C = 5 . C = 6 . C = 7 . (cid:121) (cid:88) (cid:121)(cid:121) (cid:88) (cid:56)(cid:82) (cid:88) (cid:121)(cid:82) (cid:88) (cid:56)(cid:107) (cid:88) (cid:121) (cid:121) (cid:56) (cid:82)(cid:121) (cid:82)(cid:56) (cid:107)(cid:121) (cid:105) ˜ x (cid:121)(cid:121) (cid:88) (cid:107)(cid:121) (cid:88) (cid:57)(cid:121) (cid:88) (cid:101)(cid:121) (cid:88) (cid:51)(cid:82) (cid:76)(cid:27) (cid:121) (cid:88) (cid:121)(cid:121) (cid:88) (cid:56)(cid:82) (cid:88) (cid:121)(cid:82) (cid:88) (cid:56)(cid:107) (cid:88) (cid:121) (cid:121) (cid:56) (cid:82)(cid:121) (cid:82)(cid:56) (cid:107)(cid:121) (cid:105) ˜ x (cid:121)(cid:121) (cid:88) (cid:107)(cid:121) (cid:88) (cid:57)(cid:121) (cid:88) (cid:101)(cid:121) (cid:88) (cid:51)(cid:82) (cid:76)(cid:27) (cid:121) (cid:88) (cid:121)(cid:121) (cid:88) (cid:56)(cid:82) (cid:88) (cid:121)(cid:82) (cid:88) (cid:56)(cid:107) (cid:88) (cid:121) (cid:121) (cid:56) (cid:82)(cid:121) (cid:82)(cid:56) (cid:107)(cid:121) (cid:105) ˜ x (cid:121)(cid:121) (cid:88) (cid:107)(cid:121) (cid:88) (cid:57)(cid:121) (cid:88) (cid:101)(cid:121) (cid:88) (cid:51)(cid:82) (cid:76)(cid:27) a)b)c) Figure 4: Υ for simplified model with a) C = 5 .
3, b) C = 6 .
4, c) C = 7 . τ = 54, b) τ = 89, c) τ = 143 discrete timesteps. The grey dots correspond to missingvalues which are due to the fact that we only colour the last point of the modeled subsequence. Moreover, inrare cases the auto.arima() R function is not able to find an appropriate model. Red lines mark the stableequilibria, while the dotted red line marks the unstable equilibrium.12otted red line to the unstable one. The colors ranging from dark blue to yellow represent the stability of thepoints measured by Υ and we always color the last point of the subsequence. The simulation is initializedaround the stable equilibrium x e , where short memory of the random perturbations should prevail. Asexpected, the values of Υ remains close to 0 (corresponding to a most likely AR(1) model) as long as thesimulation oscillates around the equilibrium. The timeseries eventually approaches the neighborhood of theunstable equilibrium where long memory properties are to be expected and thus higher order ARMA(p,q)models, hence larger values of Υ, are more likely. This first transition through the unstable equilibrium iswell recognized with higher values of Υ (green dots after the dotted red line).The Anderson Darling Test can be used to find the biggest τ for which we can assume that most ofthe subsequences are sampled from a normal distribution and hence trust the ARMA model selection andfitting. As shown in figure 5, for C = 6 . τ = 67 the median ofthe p-values for all subsequences is greater than the significance level 0.05. The solid line in the gray boxes τ p v a l u e
25 35 45 55 65 75 . . . . Figure 5: Boxplot of the p-values from the Anderson Darling Test for simulated time series with C = 6 . τ AD = 67 discrete timesteps.is the median of the p-values for a fixed τ while the upper and lower border of the gray boxes refer to theupper and lower quartile of the p-values. The dotted horizontal line is the significance level. We report thatthe values of τ AD given by the Anderson Darling Test are ranging from 60 to 70 discrete points for all valuesfor C .Figure 6 summarizes the Υ values obtained for different choices of τ and different values of C in abifurcation diagram. In the figure, the equilibrium solutions of the deterministic equation (3) are shown by a a) b) c) Figure 6: Bifurcation diagram of the deterministic system (red) and Υ calculated for simulated data (graydots): a) for τ = min { mean ( T i ) | i = 1 , , } − τ = τ AD (Anderson Darling Test) and c) τ = τ AD if τ AD < τ KMeans −
5, otherwise the subseries are discarded. The red full and dotted lines showrespectively the stable and unstable branches of the bifurcation diagram.red line for the considered range of values of C . This is the same diagram as shown in Fig. 1, where the upperand lower branches of the equilibrium solution correspond to the two stable equilibria, while the middle oneis the unstable equilibrium separating the two basins of attraction of the stable equilibria. A discontinuity13n the solution is visible between the upper and middle solution branches, which is due to the discontinuityintroduced by the cutoff form of the stability function. The Υ values obtained for the simulations of thestochastic system (eq. (9)) for all considered values of C are then shown as a scatter plot along with theequilibrium solution, and the darker color corresponds to higher values of Υ. As the initial condition forall simulations is taken at the lowest equilibrium values, the transitions are expected to occur between thelower and upper equilibrium branches when the system transitions from the basin of attraction of the lowestequilibrium value to that of the highest value. High values of Υ are indeed mainly found in this region ofthe diagram.Three methods are used to select the value of τ and the results associated with these window sizes areshown in figure 6. In figure 6 panel a, τ = τ KMeans − C , e.g. for C = 5 .
4, large values of Υ are occasionally inappropriately found around the upper stable branch. For small C , the potential well will be relatively steep and the system rapidly approaches the second equilibrium, sothat the detection can be too slow. Figure 6 panel b shows the results for Υ when choosing τ according tothe Anderson Darling test, denoted as τ AD . The figure is very similar to the one using τ KMeans except thatfor C ≤ . C ’s the τ ’s chosen by the Anderson Darling test are larger than the ones estimated by the K-meansalgorithm. Consequently, the local stationarity assumption may break down. For C ≥ . τ ’s given bythe AD Test are smaller than the ones of the K-means algorithm. In these cases Υ gives a good indicationfor the stability. Figure 6 panel c is a bifurcation plot showing only the time series for which τ can be chosento satisfy both the K-means and the Anderson Darling condition, i.e. τ AD < τ KMeans −
5. Here τ := τ AD is used for the analysis and timeseries that do not satisfy the condition are discarded from the analysis.In this case, we see that large values of Υ always occur between the unstable branch and the upper stablebranch, thus Υ is capable of recognizing the location of unstable equilibria for all C and stable equilibria arenever assigned a large value of Υ. Therefore, we are confident that we can apply Υ to observational data.Moreover, for this simple bistable example system, analytical results can provide the expected time takenby the system to transition from one of the local equilibria to the bifurcation point [Krumscheid et al., 2015]and can serve as a comparison to the statistical estimations of τ obtained here. As a matter of fact, for C close to the bifurcation point, the results given by the Anderson Darling test are similar to those given byanalytical calculations. In this section we apply the stability indicator Υ on observational data obtained from one site near Dumosa,Australia for which nocturnal data are selected, and from Dome C, Antarctica for which we consider thePolar winter. When we plot ∆ T over U for both sites (see figure 7) we see a clear sign of two metastablestates: one when the wind is weak and ∆ T is large and one for strong wind where ∆ T is small. The first observational dataset consists of temperature measurements from a site near Dumosa, Victoria,Australia. The site was located in a large area with mostly homogeneous and flat terrain, covered by wheatcrops, and measurements were taken during the crop season. The temperature measurements were made onthe main tower at heights of 3 and 6m and the wind measurements at 6m. The frequency of measurementsis 1 minute. Further details about the observational site can be found in Lang et al. [2018]. As we want touse data where we can expect temperature inversions to take place we exclusively use evening and nighttimedata from March until June 2013 (89 days). Each night of data results in a timeseries of 1020 discreteobservations. Similarly to the simulated data, we use the K-means algorithm and the Anderson DarlingTest to choose the window length τ . The results are shown for all nights considered together in Figure 8.14 (cid:107) (cid:57) (cid:101) (cid:51) (cid:82)(cid:121) (cid:82)(cid:107) (cid:121)(cid:82)(cid:107)(cid:106) (cid:108) (cid:101)(cid:75) (cid:104) (cid:101) (cid:75) (cid:64) (cid:104) (cid:106) (cid:75) (cid:121) (cid:107) (cid:57) (cid:101) (cid:51) (cid:82)(cid:121) (cid:82)(cid:107) (cid:121)(cid:56)(cid:82)(cid:121)(cid:82)(cid:56)(cid:107)(cid:121)(cid:107)(cid:56) (cid:108) (cid:51)(cid:75) (cid:104) (cid:78) (cid:88) (cid:57) (cid:75) (cid:64) (cid:104) (cid:98) (cid:95) (cid:89) (cid:73)(cid:51)(cid:121)(cid:113)(cid:75) (cid:64)(cid:107) (cid:95) (cid:89) (cid:61)(cid:82)(cid:121)(cid:121)(cid:113)(cid:75) (cid:64)(cid:107) a) b) Figure 7: Temperature inversion as a function of wind speed as observed at Dumosa, Australia (a) and DomeC, Antarctica (b). The color in panel b corresponds to lower and higher incoming radiation.According to these tests the maximal τ for which we can assume normality ( τ AD ) is 19 discrete observations a) b) Figure 8: a) Boxplot of the p-values from the Anderson Darling test and b) Clustered data with K-meansfor the Dumosa data, for all 89 nights.(hence 19 minutes with the sampling frequency of 1 minute) and the τ which corresponds to the minimalmean residence time in one of three clusters ( τ KMeans ) is 22. Hence we have τ AD < τ KMeans and choosing τ AD should be appropriate. The results for Υ applied to all 89 nights with both choices of window lengthsare given in figure 9. The results highlight a lower branch with low values of Υ, or dynamics identified asstable, and an upper branch with high values of Υ, or dynamics identified as unstable. In some cases, aproper ARMA(p,q) model cannot be fitted by the statistical methods, resulting in absence of results forsome windows. Generally, a reliable ARMA(p,q) fit becomes difficult for a time series with less than 20observations, and the estimated window lengths are on the lowest end to obtain the statistical estimations.Figure 10 shows the time evolution of ∆ T when conditionally averaged for all nights with the wind speed(wsp) being in a given category. The corresponding time evolution of Υ is shown for the same conditionalaverages. The window length here is chosen as the most restrictive criterion τ = τ AD < τ KMeans . Forlow wind speeds, the Υ values are high (on average), which implies that in this case we have an unstablesystem. Note that in this dataset, a leveling-off of the temperature inversion for low wind speeds (whichcould correspond to the stable equilibrium of a strong inversion according to the model of van de Wiel et al.[2017]) is not very evident from Figure 9. It could be that the temperature inversion does not have time15igure 9: Observed temperature inversion versus wind speed relation for the Dumosa data. Colored accordingto Υ with different window lengths τ : a) τ = τ AD , b) τ = τ KMeans . The grey dots correspond to missingvalues which are due to the fact that we only color the last point of the modeled subsequence. Moreover, insome cases the auto.arima() R function is not able find an appropriate model. a) b)
Figure 10: Time series of mean(∆ T ) (a) and mean(Υ) (b) for different wind speed (wsp) categories forDumosa data, calculated with τ KMeans . The shaded area is the standard deviation.to reach the stable equilibrium during the night, or that other instabilities which are not considered in thesimplified model arise in strong stability conditions. For example, flow instabilities such as submeso motions,which are favored in strongly stratified situations, could make the system dynamically unstable.
The Dome C data was measured at the Concordia Research Station which is located on the AntarcticaPlateau. It is a French-Italian research facility that was built 3233m above sea level. It is extensivelydescribed for example in Genthon et al. [2010]. The Dome C dataset contains 10-min averaged meteorologicaldata from 2017. Regimes and their transitions were analyzed by Vignon et al. [2017] and Baas et al. [2019].Important for our analysis are measurements of the temperature at height 9.4 m and surface, the wind speed(m/s) at height 8 m and the radiation made in the polar night which is from March to September. We focuson the polar night during which multiple regime transitions take place. Following van de Wiel et al. [2017]the data is classified into two subcategories of radiative forcing being the sum of net shortwave and incominglongwave radiation: R + = K ↓ − K ↑ + L ↓ . Strong cooling is favored in cases of low incoming radiation andwhen plotting ∆ T = T . m − T s over the wind speed U m a back-folding of the points becomes apparent when R + < W m − (van de Wiel et al. [2017], their Figure 6 and less clearly in our Figure 7). Therefore, we focuson the case when R + < W m − . We apply Υ to the longest consecutive time series with R + < W m − which is from 2017-08-03 10:50 to 2017-08-24 21:50, i.e. 3091 data points. Our analyzed time series is thusshorter than the one visualized in van de Wiel et al. [2017], which explains the differences in the scatter plot.16gain we choose τ with the Anderson Darling Test and the K-Means algorithm. The value for τ given bythe Anderson Darling Test τ AD = 43 observations (with a corresponding window duration of 430 minutes)is much larger than the one given by the K-means algorithm ( τ KMeans = 10 observations, corresponding toa window duration of 100 minutes), which is in fact too few points to expect a good fit for the ARMA(p,q)model. In figure 11 we see that no transitions are recognized by Υ with τ = τ AD (left) but with τ = τ KMeans (right) some transitions are noted. This indicates that the data frequency of this data set is not high enough (cid:107) (cid:57) (cid:101) (cid:51) (cid:82)(cid:121) (cid:82)(cid:107) (cid:121)(cid:56)(cid:82)(cid:121)(cid:82)(cid:56)(cid:107)(cid:121)(cid:107)(cid:56) (cid:108) (cid:51)(cid:75) (cid:104) (cid:78) (cid:88) (cid:57) (cid:75) (cid:64) (cid:104) (cid:98) (cid:107) (cid:57) (cid:101) (cid:51) (cid:82)(cid:121) (cid:82)(cid:107) (cid:121)(cid:56)(cid:82)(cid:121)(cid:82)(cid:56)(cid:107)(cid:121)(cid:107)(cid:56) (cid:108) (cid:51)(cid:75) (cid:104) (cid:78) (cid:88) (cid:57) (cid:75) (cid:64) (cid:104) (cid:98) (cid:121)(cid:121) (cid:88) (cid:56)(cid:82) (cid:76)(cid:27) a) b) Figure 11: Temperature inversion between 9.4 m and the surface, as a function of wind speed at 8m asobserved at Dome C. Colored according to Υ with different window lengths τ , expressed in number ofdiscrete observations: a) τ = 10 (K-means) and b) τ = 43 (AD Test)to give reliable results for the stability indicator. As discussed when applying the indicator Υ on the Dome C dataset, there is strong indication for the datafrequency being crucial for the reliability of the results when applying the stability indicator. Observationaldata is often stored in block averages, e.g. measurements over 5 minute time window are averaged into 1 datapoint. The issue with this can be that the data frequency can be too low to sample typical fluctuations duringthe observed transition. In more detail, if the time taken by the system to transition from one metastablestate to the next is less than approximately 20 discrete measurement points (the minimum needed to haverelevant statistical results according to our tests), then the approach may not be applicable. Therefore,data frequency needs to be high enough to give reliable results for Υ. As a comparative study to illustratethis point, we block average the temperature measurements for the Dumosa data, such that we repeat theanalysis based on 5-min averaged data instead of 1-min averages. Thereby, we reduce the length of the timeseries for each individual night from 1020 to 204 data points. Again we choose τ with the Anderson DarlingTest and the K-means algorithm. There is a clear distinction between the τ estimated by the AndersonDarling Test and the one given by the K-means algorithm for the 5 min data, contrary to the 1 min datawhere both methods suggested comparable window lengths. Indeed, τ AD = 31 and τ KMeans = 7 for the 5min averaged data whereas τ AD = 19 and τ KMeans = 22 for the 1 min data. The small value for τ givenby the K-means algorithm for the 5 min averaged data suggests that there is only a small amount of pointscovering the transition time and we cannot fit an ARMA(p,q) model properly to subsequences this short.Moreover, as the value for τ given by the Anderson Darling Test is much bigger than the amount of pointscovering the transition time we do not expect reliable results for Υ with this τ . Figure 12 confirms thishypothesis. The left plot is with τ = τ AD and the right one with τ = τ KMeans . In this study we analyzed the potential of a statistical indicator to be used to detect the system’s proximityto critical regime transitions in the near-surface temperature inversion. The statistical indicator evaluatesthe dynamical stability of time series resulting from a dynamical system and was initially suggested inFaranda et al. [2014]. Based on idealized numerical simulations, van Hooijdonk et al. [2016] had foundthe presence of early-warning signs in the turbulent flow field before a transition from weakly stable to17igure 12: Observed temperature inversion versus wind speed relation for the 5-min averaged Dumosa data.Colored according to Υ with different window lengths τ , expressed in number of discrete observations: a) τ = τ AD and b) τ = τ KMeans .strongly stable conditions. These signs included a critical slowing down, referring to the fact that dynamicalsystems tend to recover slower from perturbations when approaching a transition point in the dynamics.This slowdown was evaluated based on fluctuations of the temperature field and the early-warning signalrelied on a change in the variance. Such metrics, which are often used in studies of tipping points, canbecome problematic when the underlying dynamics is highly non-stationary, as an increase of variance couldbe due to the non-stationarity of the system without implying a transition [Lenton et al., 2012, Farandaet al., 2014]. The typical scatter of atmospheric field data and their inherent non-stationarity makes theapplication of classical critical slowdown metrics difficult. The metric presented and used here is differentin that it statistically quantifies the deviation from the dynamics expected when the system is close toa stable equilibrium. Specifically, the indicator is based on ARMA modeling with a moving window forwhich local stationarity is assumed, and the distance from stable equilibrium dynamics is evaluated basedon a Bayesian information criterion. The indicator crucially relies on an appropriate window length andwe suggested two methods to select its value in a data-driven manner. That is, both methods can be usedwhen the underlying model governing the dynamics is unknown, such that those can be applied to fielddata with significant scatter. Our suggestion to ensure reliable results is to use a combination of bothapproaches. The shortest residence time around an equilibrium estimated through the K-means approachprovides an upper bound to select a window length that respects the timescale of the system, i.e. a lengththat ensures local stationarity for ARMA model fitting. The window length should be selected as shorteror equal to this upper bound, and such that the data within individual windows mostly satisfy normality toensure reliable Bayesian inference. The Anderson Darling normality test is appropriate, but an improvementof the clustering approach to estimate the residence time around an equilibrium (here done based on asimple K-means clustering approach) would be beneficial. Based on this approach, we find that a nocturnaltemperature inversion dataset with a sampling frequency of 1 minute can be analyzed successfully using awindow length of approximately 20 minutes. Slower sampling frequency did not lead to conclusive results.The conceptual model introduced by van de Wiel et al. [2017] was developed to understand regime tran-sitions in the near-surface temperature inversion and can support scenarios with multiple stable equilibria.For our purpose of identifying the system’s proximity to regime transitions, it offers an ideal model for whichthe theoretical dynamical stability can be calculated analytically. We extended the model to include randomperturbations in the dynamics and used the resulting stochastic model to provide a test dataset on whichto evaluate the potential of the indicator of regime transitions. In this stochastic system, small-scale per-turbations can be amplified due to the nonlinearity, resulting in transitions between the bi-stable equilibria.Our simulations show such noise-induced regime transitions, successfully identified by the indicator Υ. Moreresearch would however be beneficial in order to assess the type of noise that is appropriate to representrandomized dynamics of the SBL.The application to field data was done for one nocturnal dataset and one Polar dataset. In their discus-sion, van de Wiel et al. [2017] suggest that the strength of the thermal coupling between the soil and the18tmosphere may be a key process to distinguish between cases where the temperature inversion has a uniquestable equilibria and cases with bi-stable equilibria, separated by an unstable equilibrium. The wind-speeddependence of observational scatter is partly attributed to the existence of a dynamically unstable branch inthe system in cases where the thermal coupling is weak. In both datasets considered in our analysis, a weakthermal coupling is to be expected. Clearly in the Polar dataset, the snow surface leads to a weak thermalcoupling between the atmosphere and the soil [Vignon et al., 2017, van de Wiel et al., 2017]. The nocturnaldataset originates from a wheat crop near Dumosa, Australia, probably resulting in a weak thermal couplingas well. While the Dome C data did not have the required sampling rate in order to have reliable estimatesof the dynamical stability, the Dumosa data were found to have a clear signal with one dynamically stablebranch and one dynamically unstable branch. A second dynamically stable branch corresponding to a stronginversion was not clearly observed. This data-driven result agrees with the theoretical result of van de Wielet al. [2017], namely that a dynamically unstable branch exists for a certain range of wind speeds in case ofweak atmosphere-surface thermal coupling. Note that this is an idealized model and other non-representedphysical processes may be at work and impact the interpretations. This result is nevertheless promising forthe use of the indicator as an early-warning signal of regime transitions. Extending the analysis to a Polarnight with an appropriate sampling frequency would be very interesting, as multiple regime transitions occurduring the long-lived temperature inversion [Baas et al., 2019]. Moreover, comparing results obtained for asite with strong atmosphere-surface thermal coupling would provide great insight to compare the dynamicalstability of field data to the dynamical stability predicted by the conceptual model.To be noted is the fact that the conceptual model is derived for a temperature inversion between thesurface and a height at which the wind speed stays relatively constant during the night, found to be approx-imately 40 m at Cabauw in the Netherlands and 10 m at Dome C. The Dumosa dataset did not offer thepossibility to select a height with such a constant wind speed. The measurements, taken at lower heightsin this case, will be prone to submeso-scale activity, inducing perturbations of the shallow inversion whichcould affect the dynamical stability of the time series. In fact, the earlier application of the dynamicalstability indicator to SBL data in Nevo et al. [2017] showed that higher stability corresponded to unstabledynamics of the vertical velocity fluctuations and of the wind speed. More analyses would be needed toassess the influence of the measurement height on the evaluated dynamical stability of the temperatureinversion. Nevertheless, our results encourage the use of the statistical dynamical stability as a metric todetect nearing regime transitions in the SBL. The ability to detect nearing regime transitions in atmosphericnumerical weather prediction and climate models could offer a possibility to use a different type of SBLparameterization in those specific cases without relying on the assumption of turbulence stationarity.
Acknowledgements
The authors wish to thank Christophe Genthon and Etienne Vignon for providing the Dome C data, ob-tained as part of the CALVA observation project supported by the French polar institute IPEV, and forlively discussions. The research was funded by the Deutsche Forschungsgemeinschaft (DFG) through grantnumber VE 933/2-2. A.K acknowledges funding from the DFG through the Collaborative Research CenterCRC1114, ”Scaling Cascades in Complex system”. The scientific exchanges between N.V and D.F weregreatly facilitated by funding through the DAAD exchange program Procope through the project ”Data-driven dynamical stability of stably stratified turbulent flows”. We are grateful to anonymous reviewers forinsightful comments that helped us improve the manuscript.
Appendix
The clustering is done using the K-means algorithm with the following steps: • Input: k = x i − τ +1 , . . . , x i • Place centroids c , . . . , c k at random locations. 19 Repeat until none of the cluster assignments change: – for each point x i find nearest centroid c j and assign x i to cluster j – for each cluster j = 1 , . . . , k calculate new centroid c j = mean of all points x i assigned to cluster j in previous step. The the Anderson-Darling (AD) normality test statistic is based on the squared difference between theempirical distribution function estimated based on the sample, F n ( x ), and the normal distribution F ∗ ( x ).The statistic for this test is, W n = n (cid:90) ∞−∞ [ F n ( x ) − F ∗ ( x )] ψ ( F ∗ ( x )) dF ∗ ( x )where ψ is a non-negative weight function which is used to emphasize the tails of the presumed distribution.We use the modified AD statistic given by D’Agostino and Stephens (1986) which takes into accounts thesample size n W ∗ n = W n (1 + 0 . /n + 2 . /n ) . The null hypothesis of this normality test is that the data are sampled from a normal distribution. Whenthe p-value is greater than the predetermined critical value ( α = 0 . The full procedure to apply the statistical indicator to a timeseries follows the following steps: • Input: Timeseries x , . . . , x T • Evaluate the window length τ Kmeans : – Cluster the timeseries in k clusters (expected number of equilibrium states) using the K-meansalgorithm. – For each cluster, calculate the mean residence time of the timeseries. – The minimal mean residence time over the k clusters provides an upper bound for τ Kmeans . • Evaluate the window length τ AD : – For a range of window lengths τ min < τ < τ max , apply the AD test statistic to the timeseries ina moving window approach. – For each τ , calculate the median of the p-values obtained over all windows. – Select the largest τ for which the median p-value is greater than the predetermined critical value( α = 0 .
05) as τ AD . • If τ AD < τ Kmeans , select τ AD as a window length. Else, ARMA model fitting may be innappropriate. • Repeat for each window of length τ AD : – Select the best fitting ARMA(p,q) model (minimal BIC). – Fit an ARMA(1,0) to the window and calculate the BIC. – Calculate Υ. High values will indicate transitions. The values themselves may depend on thedataset. 20 eferences
Carsten Abraham and Adam H Monahan. Climatological Features of the Weakly and Very Stably StratifiedNocturnal Boundary Layers. Part I: State Variables Containing Information about Regime Occupation.
Journal of Atmospheric Sciences , 76(11):3455–3484, November 2019.Ot´avio C Acevedo, Osvaldo L L Moraes, Gerv´asio A Degrazia, and Luiz E Medeiros. Intermittency and theExchange of Scalars in the Nocturnal Surface Layer.
Boundary-Layer Meteorology , 119(1):41–55, March2006.Ot´avio C Acevedo, Larry Mahrt, Franciano S Puhales, Felipe D Costa, Luiz E Medeiros, and Gerv´asio ADegrazia. Contrasting structures between the decoupled and coupled states of the stable boundary layer.
Quarterly Journal of the Royal Meteorological Society , 142(695):693–702, December 2015.Ot´avio C Acevedo, Rafael Maroneze, Felipe D Costa, Franciano S Puhales, Gerv´asio A Degrazia, Luis G NMartins, Pablo E S Oliveira, and Luca Mortarini. The Nocturnal Boundary Layer Transition from Weaklyto Very Stable. Part 1: Observations.
Quarterly Journal of the Royal Meteorological Society , 145(725):3577–3592, August 2019.T. W. Anderson and D. A. Darling. Asymptotic theory of certain ”goodness of fit” criteria based on stochasticprocesses.
The Annals of Mathematical Statistics , 23(2):193–212, 1952.Peter Baas, Bas J H van de Wiel, S J A Linden, and Fred C Bosveld. From Near-Neutral to StronglyStratified: Adequately Modelling the Clear-Sky Nocturnal Boundary Layer at Cabauw.
Boundary-LayerMeteorology , 166(2):217–238, October 2017.Peter Baas, Bas J H van de Wiel, Erik van Meijgaard, Etienne Vignon, Christophe Genthon, Steven J Avan der Linden, and Stephan R de Roode. Transitions in the wintertime near-surface temperature inversionat Dome C, Antarctica.
Quarterly Journal of the Royal Meteorological Society , 145(720):930–946, March2019.P. J. Brockwell and R. A. Davis.
Introduction to Time Series and Forecasting . Springer, 3. edition, 2016.ISBN 978-3-319-29852-8.Luonan Chen, Rui Liu, Zhi-Ping Liu, Meiyi Li, and Kazuyuki Aihara. Detecting early-warning signals forsudden deterioration of complex diseases by dynamical network biomarkers.
Scientific reports , 2:342, 2012.Judith M M Donda, Ivo G S van Hooijdonk, Arnold F Moene, Harmen J J Jonker, G J F van Heijst, HermanJ H Clercx, and Bas J H van de Wiel. Collapse of turbulence in stably stratified channel flow: a transientphenomenon.
Quarterly Journal of the Royal Meteorological Society , 141(691):2137–2147, July 2015.D. Faranda, B. Dubrulle, and F. M. E. Pons. Statistical early-warning indicators based on autoregressivemoving-average models.
Journal of Physics A: Mathematical and Theoretical , 47(25):252001, Jun 2014.URL https://doi.org/10.1088/1751-8113/47/25/252001 .C. Genthon, M. S. Town, D. Six, V. Favier, S. Argentini, and A. Pellegrini. Meteorological atmosphericboundary layer measurements and ecmwf analyses during summer at dome c, antarctica.
Journal ofGeophysical Research: Atmospheres , 115(D5), 2010. doi: 10.1029/2009JD012741.J. A. Hartigan and M. A. Wong. Algorithm as 136: A k-means clustering algorithm.
Journal of the RoyalStatistical Society. Series C (Applied Statistics) , 28(1):100–108, 1979.Alan Hastings and Derin B Wysham. Regime shifts in ecological systems can occur with no warning.
Ecologyletters , 13(4):464–472, 2010.Amber M Holdsworth and Adam H Monahan. Turbulent Collapse and Recovery in the Stable BoundaryLayer Using an Idealized Model of Pressure-Driven Flow with a Surface Energy Budget.
Journal ofAtmospheric Sciences , 76(5):1307–1327, May 2019.21AM Holtslag, M Tjernstr¨om, S Basu, A C M Beljaars, Gunilla Svensson, P Baas, B Beare, Fred C Bosveld,Joan Cuxart, J Lindvall, G J Steeneveld, and Bas J H van de Wiel. Stable Atmospheric Boundary Layersand Diurnal Cycles: Challenges for Weather and Climate Models.
Bulletin of the American MeteorologicalSociety , 94(11):1691–1706, November 2013.R. Hyndman, G. Athanasopoulos, C. Bergmeir, and et. al. Forecasting functions for time series and linearmodels, 2019. URL http://pkg.robjhyndman.com/forecast . R package version 8.4.S. Krumscheid, M. Pradas, G. A. Pavliotis, and S. Kalliadasis. Data-driven coarse graining in ac-tion: Modeling and prediction of complex systems.
Phys. Rev. E , 92:042139, Oct 2015. URL https://link.aps.org/doi/10.1103/PhysRevE.92.042139 .Yuri A Kuznetsov.
Elements of applied bifurcation theory , volume 112. Springer Science & Business Media,2013.Francisco Lang, Danijel Beluˇsi´c, and Steven Siems. Observations of Wind-Direction Variability in the Noc-turnal Boundary Layer.
Boundary-Layer Meteorology , 166(1):51–68, September 2018.Margaret A LeMone, Wayne M Angevine, Christopher S Bretherton, Fei Chen, Jimy Dudhia, Evgeni Fe-dorovich, Kristina B Katsaros, D H Lenschow, Larry Mahrt, EG Patton, Jielun Sun, Michael Tjernstr¨om,and Jeffrey Weil. 100 Years of Progress in Boundary Layer Meteorology.
Meteorological Monographs , 59:9.1–9.85, January 2018.TM Lenton, VN Livina, V Dakos, EH Van Nes, and M Scheffer. Early warning of climate tipping pointsfrom critical slowing down: comparing methods to improve robustness.
Philosophical Transactions of theRoyal Society A: Mathematical, Physical and Engineering Sciences , 370(1962):1185–1204, 2012.Larry Mahrt. Stably Stratified Atmospheric Boundary Layers.
Annual Review of Fluid Mechanics , 46:23–45,2014.Rafael Maroneze, Ot´avio C Acevedo, Felipe D Costa, Franciano S Puhales, Giuliano Demarco, and LucaMortarini. The Nocturnal Boundary Layer Transition from Weakly to Very Stable. Part 2: NumericalSimulation with a Second Order Model.
Quarterly Journal of the Royal Meteorological Society , 145(725):3593–3608, August 2019.R T McNider. Predictability of the stable atmospheric boundary layer.
Journal of the atmospheric sciences ,52(10):1602–1614, May 1995.Adam H Monahan, Tim Rees, Yanping He, and Norman McFarlane. Multiple Regimes of Wind, Stratifica-tion, and Turbulence in the Stable Boundary Layer.
Journal of Atmospheric Sciences , 72(8):3178–3198,August 2015.G. Nevo, N. Vercauteren, A. Kaiser, B. Dubrulle, and D. Faranda. Statistical-mechanical approach to studythe hydrodynamic stability of the stably stratified atmospheric boundary layer.
Phys. Rev. Fluids , 2(8):084603, Aug 2017. doi: 10.1103/PhysRevFluids.2.084603.Irina Sandu, A C M Beljaars, Peter Bechtold, Thorsten Mauritsen, and Gianpaolo Balsamo. Why is it sodifficult to represent stably stratified conditions in numerical weather prediction (NWP) models?
Journalof Advances in Modeling Earth Systems , 5(2):117–133, June 2013.Marten Scheffer, Jordi Bascompte, William A Brock, Victor Brovkin, Stephen R Carpenter, Vasilis Dakos,Hermann Held, Egbert H Van Nes, Max Rietkerk, and George Sugihara. Early-warning signals for criticaltransitions.
Nature , 461(7260):53, 2009.Marten Scheffer, Stephen R Carpenter, Timothy M Lenton, Jordi Bascompte, William Brock, Vasilis Dakos,Johan Van de Koppel, Ingrid A Van de Leemput, Simon A Levin, Egbert H Van Nes, et al. Anticipatingcritical transitions. science , 338(6105):344–348, 2012.G. Schwarz. Estimating the dimension of a model.
The Annals of Statistics , 6(2):461–464, 1978.22. A. Stephens. Edf statistics for goodness of fit and some comparisons.
Journal of the Amer-ican Statistical Association , 69(347):730–737, 1974. doi: 10.1080/01621459.1974.10480196. URL .Jielun Sun, Larry Mahrt, Robert M Banta, and Yelena L Pichugina. Turbulence Regimes and TurbulenceIntermittency in the Stable Boundary Layer during CASES-99.
Journal of Atmospheric Sciences , 69(1):338–351, January 2012.Jielun Sun, D H Lenschow, Margaret A LeMone, and Larry Mahrt. The Role of Large-Coherent-Eddy Trans-port in the Atmospheric Surface Layer Based on CASES-99 Observations.
Boundary-Layer Meteorology ,160(1):83–111, July 2016.D J Thomson. Criteria for the selection of stochastic models of particle trajectories in turbulent flows.
Journal of Fluid Mechanics , 180(-1):529–556, July 1987.Bas J H van de Wiel, Arnold F Moene, G J Steeneveld, O K Hartogensis, and AAM Holtslag. Predictingthe Collapse of Turbulence in Stably Stratified Boundary Layers.
Flow, Turbulence and Combustion , 79(3):251–274, 2007.Bas J H van de Wiel, A F Moene, Harmen J J Jonker, P Baas, S Basu, Judith M M Donda, Jielun Sun, andAAM Holtslag. The Minimum Wind Speed for Sustainable Turbulence in the Nocturnal Boundary Layer.
Journal of Atmospheric Sciences , 69(11):3116–3127, November 2012a.Bas J H van de Wiel, Arnold F Moene, and Harmen J J Jonker. The Cessation of Continuous Turbulenceas Precursor of the Very Stable Nocturnal Boundary Layer.
Journal of Atmospheric Sciences , 69(11):3097–3115, November 2012b.Bas J H van de Wiel, Etienne Vignon, Peter Baas, Ivo G S van Hooijdonk, Steven J A van der Linden,J Antoon van Hooft, Fred C Bosveld, Stefan R de Roode, Arnold F Moene, and Christophe Genthon.Regime Transitions in Near-Surface Temperature Inversions: A Conceptual Model.
Journal of AtmosphericSciences , 74(4):1057–1073, April 2017.Ivo G S van Hooijdonk, Judith M M Donda, Herman J H Clercx, Fred C Bosveld, and Bas J H van de Wiel.Shear Capacity as Prognostic for Nocturnal Boundary Layer Regimes.
Journal of Atmospheric Sciences ,72(4):1518–1532, April 2015.Ivo G S van Hooijdonk, M Scheffer, H J H Clercx, and Bas J H van de Wiel. Early Warning Signals forRegime Transition in the Stable Boundary Layer: A Model Study.
Boundary-Layer Meteorology , 162(2):1–24, October 2016.Nikki Vercauteren and Rupert Klein. A Clustering Method to Characterize Intermittent Bursts of Turbulenceand Interaction with Submesomotions in the Stable Boundary Layer.
Journal of the atmospheric sciences ,72(4):1504–1517, 2015.Nikki Vercauteren, Larry Mahrt, and Rupert Klein. Investigation of interactions between scales of motionin the stable boundary layer.
Quarterly Journal of the Royal Meteorological Society , 142(699):2424–2433,jun 2016.Nikki Vercauteren, Vyacheslav Boyko, Amandine Kaiser, and D Beluˇsi´c. Statistical Investigation of FlowStructures in Different Regimes of the Stable Boundary Layer.
Boundary-Layer Meteorology , 173(2):143–164, November 2019.Etienne Vignon, Bas J H van de Wiel, Ivo G S van Hooijdonk, Christophe Genthon, Steven J A van derLinden, J Antoon van Hooft, Peter Baas, William Maurel, Olivier Traull´e, and Giampietro Casasanta.Stable boundary-layer regimes at Dome C, Antarctica: observation and analysis.