Science, Faculty of
Mathematics, Department of
DSpace
UBCV
Yin, Ziming
Master of Science - MSc
University of British Columbia
The Atlantic Ocean circulation transports the warmth and salinity northwards and it is a crucial part of the world climate. The flipping of this circulation is hypothesized to be the explanation of several abrupt climate changes over the past 2.6 million years. The pattern of this circulation can be determined from ocean models, and a typical pattern, which is also used in this thesis, is described by the Stommel model. There are two bifurcations in the Stommel model — one is a smooth bifurcation and the other is a non-smooth bifurcation. When the system approaches either bifurcation, the presence of noise makes the system escape from the current state easily. However, the properties of these two bifurcations are widely different. Here for the smooth bifurcation, the probability of tipping is judged by some well-known linear and nonlinear methods, including the linear decay rate, the potential energy and the skewness. Whereas for the non-smooth bifurcation, if there is no discernable trend from the above methods, we find out that one is also able to obtain the early warning signal by a unique feature of the Stommel model: the eigenvalues of the system matrix when approaching the non-smooth bifurcation are complex numbers, which lead to periodic circulations. The change of the periodic circulation can be observed by the variation of its power spectral density, which is important in showing the proximity to tipping. In this thesis, the coherence resonance is used to define this variation. After some extensive studies, we find that the combination of the above methods are able to give us an early warning signal as the system moves closer to the bifurcations.
Noise Induced Tipping in Ocean Circulation Based on theStommel ModelbyZiming YinB.Sc., Beijing Jiaotong University, 2016A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFMaster of ScienceinTHE FACULTY OF GRADUATE AND POSTDOCTORALSTUDIES(Mathematics)The University of British Columbia(Vancouver)August 2018c© Ziming Yin, 2018The following individuals certify that they have read, and recommend to the Fac-ulty of Graduate and Postdoctoral Studies for acceptance, a thesis entitled:Noise Induced Tipping in Ocean Circulation Based on the Stommel Modelsubmitted by Ziming Yin in partial fulfillment of the requirements of the degree ofMaster of Science in Applied Mathematics.Examining Committee:Rachel Kuske, Professor of MathematicsSupervisorWayne NagataSupervisory Committee MemberiiAbstractThe Atlantic Ocean circulation transports the warmth and salinity northwards andit is a crucial part of the world climate. The flipping of this circulation is hypothe-sized to be the explanation of several abrupt climate changes over the past 2.6 mil-lion years. The pattern of this circulation can be determined from ocean models,and a typical pattern, which is also used in this thesis, is described by the Stommelmodel. There are two bifurcations in the Stommel model — one is a smooth bi-furcation and the other is a non-smooth bifurcation. When the system approacheseither bifurcation, the presence of noise makes the system escape from the currentstate easily. However, the properties of these two bifurcations are widely differ-ent. Here for the smooth bifurcation, the probability of tipping is judged by somewell-known linear and nonlinear methods, including the linear decay rate, the po-tential energy and the skewness. Whereas for the non-smooth bifurcation, if thereis no discernable trend from the above methods, we find out that one is also ableto obtain the early warning signal by a unique feature of the Stommel model: theeigenvalues of the system matrix when approaching the non-smooth bifurcation arecomplex numbers, which lead to periodic circulations. The change of the periodiccirculation can be observed by the variation of its power spectral density, which isimportant in showing the proximity to tipping. In this thesis, the coherence reso-nance is used to define this variation. After some extensive studies, we find thatthe combination of the above methods are able to give us an early warning signalas the system moves closer to the bifurcations.iiiLay SummaryWe study the overturning behavior inside the Atlantic Ocean circulation, whichis responsible for some of the abrupt climate shifts. Inside the circulation, thereare a smooth and a non-smooth changes, which are also called bifurcations in thisthesis. Our main task is to predict the smooth and non-smooth bifurcations ahead,so that it can help people understand the overturning behavior inside the AtlanticOcean circulation. Therefore we test some well-established methods on the twobifurcations. Because the result of the non-smooth bifurcation is not solid, we alsoexplore some other schemes to make up for the insufficient warning information.ivPrefaceThis thesis is original, unpublished, independent work by the author, Ziming Yin.vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xviii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Bifurcation, tipping and hysteresis . . . . . . . . . . . . . . . . . 11.2 The Stommel model . . . . . . . . . . . . . . . . . . . . . . . . . 41.3 Application—tipping in the Atlantic ocean . . . . . . . . . . . . . 72 Background on the Stommel Model . . . . . . . . . . . . . . . . . . 92.1 The static case . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92.1.1 Bifurcations in the Stommel model . . . . . . . . . . . . 92.1.2 Location of the smooth saddle node bifurcation . . . . . . 112.1.3 Behavior near the equilibrium . . . . . . . . . . . . . . . 112.2 Varying parameter . . . . . . . . . . . . . . . . . . . . . . . . . . 152.3 Possible noise sources . . . . . . . . . . . . . . . . . . . . . . . . 16vi2.4 Simplified 1D Hysteresis model . . . . . . . . . . . . . . . . . . 182.4.1 Bifurcations in the 1D model . . . . . . . . . . . . . . . 193 Noise-induced tipping near a smooth saddle node bifurcation . . . . 213.1 Time series analysis methods . . . . . . . . . . . . . . . . . . . . 223.1.1 Gaussian kernel detrending . . . . . . . . . . . . . . . . . 223.1.2 Estimate of the linear decay rate . . . . . . . . . . . . . . 243.1.3 Nonlinearity—information from the potential well . . . . 313.1.4 Skewness . . . . . . . . . . . . . . . . . . . . . . . . . . 374 1-D model with smooth and non-smooth bifurcations . . . . . . . . . 394.1 Location of the smooth saddle-node bifurcation . . . . . . . . . . 414.1.1 Time series and Gaussian kernel detrending . . . . . . . . 414.1.2 Estimate the linear decay rate . . . . . . . . . . . . . . . 424.1.3 Nonlinearity information from the potential well . . . . . 454.1.4 Skewness . . . . . . . . . . . . . . . . . . . . . . . . . . 474.2 Transition for the non-smooth bifurcation . . . . . . . . . . . . . 484.2.1 Time series and Gaussian kernel detrending . . . . . . . . 484.2.2 Linear decay rate . . . . . . . . . . . . . . . . . . . . . . 494.2.3 Nonlinearity information from the potential well . . . . . 514.2.4 Skewness . . . . . . . . . . . . . . . . . . . . . . . . . . 564.3 Comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 585 The 2D Stommel model with smooth and non-smooth bifurcations . 615.1 Location of the smooth saddle node bifurcation . . . . . . . . . . 625.1.1 Time series and Gaussian kernel detrending . . . . . . . . 625.1.2 Linear Decay Rate . . . . . . . . . . . . . . . . . . . . . 635.1.3 Nonlinearity information from the potential well . . . . . 665.1.4 Skewness . . . . . . . . . . . . . . . . . . . . . . . . . . 685.2 Location of the non-smooth bifurcation . . . . . . . . . . . . . . 695.2.1 Time series and Gaussian kernel detrending . . . . . . . . 695.2.2 Estimate the linear decay rate . . . . . . . . . . . . . . . 705.2.3 Nonlinearity information from the potential well . . . . . 725.2.4 Skewness . . . . . . . . . . . . . . . . . . . . . . . . . . 78vii6 Prediction of non-smooth bifurcation via coherence resonance . . . 816.1 The coherence resonance in different realizations . . . . . . . . . 836.1.1 η1 = 3,η3 = 0.3 . . . . . . . . . . . . . . . . . . . . . . 836.1.2 η1 = 6,η3 = 0.4 . . . . . . . . . . . . . . . . . . . . . . 876.1.3 η1 = 8,η3 = 0.3 . . . . . . . . . . . . . . . . . . . . . . 896.2 Coherence resonance for smaller noise . . . . . . . . . . . . . . . 917 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 947.1 Summary of results . . . . . . . . . . . . . . . . . . . . . . . . . 947.2 Future directions . . . . . . . . . . . . . . . . . . . . . . . . . . 96Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98A Supporting Materials . . . . . . . . . . . . . . . . . . . . . . . . . . 102A.1 Noise amplitude and drift rate chosen analysis for different models 102viiiList of TablesTable 2.1 Parameters of the box model . . . . . . . . . . . . . . . . . . 17ixList of FiguresFigure 1.1 An example of bifurcation diagram generated by the static pa-rameter η2. η2 is the input parameter and Ψ is the system re-sponse. The solid line is the stable steady state, and the dottedline is the unstable steady state. See the Fig. 2.1 for detail. . . 2Figure 1.2 Two examples of tipping. The left figure is generated withoutnoise and the right figure is generated with noise. η2 is theinput parameter and Ψ is the system response. The black dot-ted line is the equilibrium of the underlying system (same asFig. 1.1). The red line describes the system when η2 is slowlyincreasing. See the Fig. 2.6 for detail. . . . . . . . . . . . . . 3Figure 1.3 One example of hysteresis behavior. Left: increasing η2. Right:decreasing η2. The red line describes how Ψ changes with η2.See the Fig. 2.5 for detail. . . . . . . . . . . . . . . . . . . . 4Figure 1.4 two-box Stommel model . . . . . . . . . . . . . . . . . . . . 5Figure 1.5 The bifurcation diagram of the Stommel model when η1 = 3and η3 = 0.3. The solid blue line is the stable solution and thedotted black line is the unstable solution. L1 marks the pointof the static non-smooth bifurcation and L2 marks the staticsmooth bifurcation. . . . . . . . . . . . . . . . . . . . . . . . 8Figure 2.1 The bifurcation diagram of the Stommel model when η1 = 3and η3 = 0.3. The solid blue line is the stable solution and thedotted black line is the unstable solution. . . . . . . . . . . . 10xFigure 2.2 The bifurcation diagram of the Stommel model when η1 = 3and η3 = 0.1,0.3,0.5,0.8 respectively. . . . . . . . . . . . . . 13Figure 2.3 The steady state of Eq (2.2). The dotted black state is unstable,and the red state is stable. η1 = 3, η3 = 0.3. . . . . . . . . . . 13Figure 2.4 The steady state of Eq (2.3). The red solid line is the stablesolution. η1 = 3, η3 = 0.3. . . . . . . . . . . . . . . . . . . . 15Figure 2.5 The bifurcation diagram when η2 is slowing varying. Left: η2increases. Right: η2 decreases. η1 = 3, η3 = 0.3. . . . . . . . 16Figure 2.6 One realization generated from Eq (2.6). η2 is slowing increas-ing. η1 = 3, η3 = 0.3, dt = 0.1, ε = 0.01, σ1 = σ2 = 0.1. . . . 18Figure 2.7 The bifurcation diagram of the 1-D hysteresis model (2.8). Thesolid line is the stable steady state and the dotted line is theunstable steady state. . . . . . . . . . . . . . . . . . . . . . . 19Figure 2.8 The steady state of process x˙ = (−µ+1)− (x−1)2. . . . . . 19Figure 2.9 The steady state of process x˙ = (−µ−1)+(x−1)2 . . . . . . 20Figure 3.1 Bifurcation diagram of the saddle node form x˙ = µ+ x2. . . . 21Figure 3.2 In both plots, the black curve represents the time series, thered line is the theoretical equilibrium, and the blue curve is theresult of filtering the time series with a Gaussian kernel. Theleft plot is obtained from the process (3.1), and the right plotis from the process (3.2). Parameters: µ(0) = 2, ∆t = 0.05,σ = 0.35, ε = 0.009, bandwidth d = 8. . . . . . . . . . . . . 23Figure 3.3 The linear decay rate for Eq (3.1) with different window lengths.The black line is the kACF , the red line represents kU as definedbelow. We try four different window lengths, and N is thetotal number of our data. Top left: the decay rate with win-dow length N/2, top right: N/3, bottom left: N/4, bottomright: N/5. Other parameters: µ(0) = 2, ∆t = 0.05, σ = 0.35,ε = 0.009, bandwidth d = 8. . . . . . . . . . . . . . . . . . . 27xiFigure 3.4 The right one is the result of the linear decay rate comingfrom the process (3.1) and the left one is coming from the pro-cess (3.2). The black curve represents kU and the red curve iskACF , the smooth blue line is the theoretical decay rate exp(−2√µ∆t)calculated from the linearized deterministic equation. The win-dow length is w = N/3. . . . . . . . . . . . . . . . . . . . . . 31Figure 3.5 The left figure is the potential well of the time series (3.1). Theright one is the linear time series (3.2). The grey surface is thepotential well estimated by Eq (3.21), the red line is a parabolaapproximated to the potential well via Eq (3.24), and the bluecurve is the time series. . . . . . . . . . . . . . . . . . . . . . 36Figure 3.6 For nonlinear process (3.1): the black dotted line is the em-pirical potential well, the red solid line is the approximatedparabola to the empirical potential well, and the grey solid lineis the theoretical potential energy. . . . . . . . . . . . . . . . 37Figure 3.7 The upper plot is the skewness γ for the process (3.1) and thelower plot is for the process (3.2). Skewness is calculated from100 realizations. . . . . . . . . . . . . . . . . . . . . . . . . 38Figure 4.1 The bifurcation diagram of the 1-D model. The dotted line isthe unstable steady state. The solid line is the stable steadystate. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39Figure 4.2 The hysteresis behavior in the 1-D model. The bifurcation di-agram when υ is slowly varying. Left: when υ is increasing.Right: when υ is decreasing. . . . . . . . . . . . . . . . . . 40Figure 4.3 The time series coming from the process (4.1) where υ is slowlyincreasing. The red line is the equilibrium, the green line is theresult of filtering the time series with a Gaussian kernel, andthe black line is the system response ys. Parameters: σ = 0.35,δ = 0.009, ∆t = 0.05, bandwidth d = 8. . . . . . . . . . . . . 41xiiFigure 4.4 The detrended data y˜s of the time series (black) in Fig. 4.3,obtained by subtracting the moving average y˜s (green line inFig 4.3) from the original time series ys (black line) beforetipping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42Figure 4.5 The linear decay rate of the detrended time series y˜s (in Fig. 4.4)with window length w=N/4: the red line is λ sACF and the blueline is λ sU . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45Figure 4.6 The nonlinear potential energy surfaces of the time series inFig. 4.3 using sliding window w=N/4. Us =−σ2emplog(pemp/2)where σemp is estimated from Equation (3.22). The grey sur-face is the potential energy, the red line is a parabola approx-imated to the potential well via Eq (3.24), the blue shadowbelow is the time series. . . . . . . . . . . . . . . . . . . . . 47Figure 4.7 The estimated skewness of process (4.1) using 100 realizationsthat run from υ = −1 to υ = 0.9. Parameters: σ = 0.35, δ =0.009, ∆t = 0.05. . . . . . . . . . . . . . . . . . . . . . . . . 48Figure 4.8 [Left] The time series coming from the 1-D model Eq (4.1)with additive noise. [Right] The detrended data y˜ns over Gaus-sian Kernel. Parameters: υ(0) = 2, σ = 0.35, δ = −0.009,∆t = 0.05, bandwidth d = 8. . . . . . . . . . . . . . . . . . . 49Figure 4.9 The linear decay rate of the detrended time series in Fig. 4.8:λ nsACF (red) from the autocorrelation function and λnsU (blue)from the potential well. The window length for the linear andnonlinear analysis is w = N/4. . . . . . . . . . . . . . . . . 51Figure 4.10 The potential well of the time series in Fig. 4.8. The greysurface is the potential well, the red line is the approximatedparabolic fit to the potential well in the last time window, andthe blue shadow below them represents the time series yns. Thepotential energy Uns = −σ2emplog(pemp/2) where σemp is esti-mated from linear decay rate λ nsU and λ nsACF . . . . . . . . . . . 53xiiiFigure 4.11 The five potential curves of time series in Fig. 4.8 (left) and inFig. 4.3 (right). Red line—the theoretical potential, Eq (4.17)for the left one and Eq (4.18) for the right one. Black line—theempirical potential curve, grey curve—the original time seriesgenerated from Eq (4.1), blue line—the theoretical equilibriumof the time series. . . . . . . . . . . . . . . . . . . . . . . . . 54Figure 4.12 The theoretical and empirical potential well: the left is for thetime series y˙ns = −υ + 2|yns|− yns|yns|+σWt and the right isfor y˙=−υ−2y+y2+σWt . Red line—the theoretical potentialcurve, black line—the empirical potential curve, grey curve—the original time series, blue line—the theoretical equilibrium. 56Figure 4.13 The estimated skewness of yns using 100 realizations generatedby process (4.1) from υ = 2 to υ = 0.5. Parameters: υ(0) = 2,σ = 0.35, δ =−0.009, ∆t = 0.05. . . . . . . . . . . . . . . . 57Figure 4.14 1st row: the left figure is the time series y˙ns = −υ − 2|yns| −yns|yns|+σWt and the right figure is y˙ =−υ−2y+ y2+σWt .2nd row: the skewness of above process. 3rd row: the proba-bility density of the process generated at fixed υ : 0.4 and 0.5. 58Figure 5.1 The two time series coming from the process (5.1) with slowlyincreasing η2. For the left time series, the noise amplitudesσ1 = σ2 = 0.2. For the right time series, σ1 = σ2 = 0.04. Theblack curve is the time series, the red line is the theoreticalequilibrium and the green line is the result of filtering the timeseries with a Gaussian kernel. Parameters: η2(0) = 0, η1 = 3,η3 = 0.3, ε = 0.005, bandwidth d = 8. . . . . . . . . . . . . . 62Figure 5.2 The linear decay rate of process (5.1) with larger noise—(a),and smaller noise—(b). Window length w = N/3 in the leftone and w = N/10 in the right one. . . . . . . . . . . . . . . 66xivFigure 5.3 The potential wells of the process (5.1). The left figure corre-sponds to the left time series in Fig. 5.1, and the right figurecorresponds to the right time series in Fig. 5.1. Here Uemp =−σ2ψ log(pemp/2)where σψ is estimated based on Equation (3.22).Sliding window w = N/3 in the left figure, and w = N/10 inthe right figure. . . . . . . . . . . . . . . . . . . . . . . . . . 67Figure 5.4 The skewness of the process (5.1) with slowly increasing η2.The left figure is generated with the larger noise and the rightfigure is generated with the smaller noise. . . . . . . . . . . . 68Figure 5.5 The two time series coming from the process (5.1) when η2 isslowly decreasing. The left figure has the noise amplitude σ1 =σ2 = 0.2, and the right figure has σ1 = σ2 = 0.04. Parameters:η2(0) = 2, η1 = 3, η3 = 0.3, ε = 0.005, ∆t = 0.05. . . . . . . 69Figure 5.6 The linear decay rate of the detrended time series ψ˜ns. (a) isfrom the time series with the larger noise and (b) is from thesmaller noise. . . . . . . . . . . . . . . . . . . . . . . . . . 72Figure 5.7 The potential well of the process (5.1). The grey surface is thepotential well, the red line is the approximated parabolic fit tothe potential well in the last time window, and the blue shadowbelow them represents the time series ψns1 and ψns2 that showin Fig. 5.5. . . . . . . . . . . . . . . . . . . . . . . . . . . . 73Figure 5.8 The theoretical (red) potential well and the empirical (black)potential well. . . . . . . . . . . . . . . . . . . . . . . . . . . 75Figure 5.9 The theoretical potential well of ψns of the deterministic ver-sion of Eq (5.1) (red) and Eq (5.9)(blue) . . . . . . . . . . . . 76Figure 5.10 The left panel is the theoretical potential well of Tns. The redline is when ψns < 0, and the blue line is when ψns > 0. Theright panel is the empirical potential well of Tns (grey) and theapproximating parabola (red). . . . . . . . . . . . . . . . . . 78Figure 5.11 The skewness of ψns of the process (5.1) with slowly decreas-ing η2. The top panel is generated with the larger noise and thebottom panel is generated with the smaller noise. . . . . . . . 79Figure 5.12 Comparison between the skewness of ψns1 and Tns1. . . . . . . 80xvFigure 5.13 Trajectories in the direction of ψ and T on the lower branch atη2 = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80Figure 6.1 The oscillation of ψns generated with fixed η2 and its corre-sponding power spectral density. The figure above is gener-ated with η2 = 1.3, which has eigenvalues −0.9± 1.3i , andthe figure below is generated with η2 = 5, which has eigenval-ues−3.0±0.8i. Other parameters: η1 = 3, η3 = 0.3, σ = 0.05,∆t = 0.05, tend = 30. . . . . . . . . . . . . . . . . . . . . . . 82Figure 6.2 (a) The bifurcation diagram of ψ when η1 = 3 and η3 = 0.3.(b) The measure of coherence resonance β of the stochasticdynamics (1.9). Other parameters: tend = 50, ∆t = 0.01, noiseamplitude σ = 0.05. . . . . . . . . . . . . . . . . . . . . . . 84Figure 6.3 The probability densities of ψns computed using 500 realiza-tions at (a) η2 = 1.2, (b) η2 = 1.067, (c) η2 = 1. Other param-eters: η1 = 3, η3 = 0.3. . . . . . . . . . . . . . . . . . . . . . 85Figure 6.4 The power spectral density of ψns at (a) η2 = 1.2, (b) η2 =1.067, (c) η2 = 1. . . . . . . . . . . . . . . . . . . . . . . . . 86Figure 6.5 (a) The bifurcation diagram when η1 = 6 and η3 = 0.4. (b)The measure of coherence resonance β computed by Eq (6.1).Other parameters: tend = 50, ∆t = 0.01, noise amplitude σ =0.05. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87Figure 6.6 The probability densities of ψns computed using 500 realiza-tions at (a) η2 = 3, (b) η2 = 2.72, (c) η2 = 2.55. Other param-eters: η1 = 6, η3 = 0.4, tend = 50. . . . . . . . . . . . . . . . 87Figure 6.7 The power spectral density of ψns at (a) η2 = 3, (b) η2 = 2.72,(c) η2 = 2.55. . . . . . . . . . . . . . . . . . . . . . . . . . . 88Figure 6.8 (a) The bifurcation diagram when η1 = 8 and η3 = 0.3. (b) Themeasure of coherence resonance β of ψns of the process (1.9)computed by Eq (6.1). Other parameters: tend = 50, ∆t = 0.01,noise amplitude σ = 0.05. . . . . . . . . . . . . . . . . . . . 89xviFigure 6.9 The probability densities (PDF) of ψns computed using 500realizations at (a) η2 = 3, (b) η2 = 2.84, (c) η2 = 2.7. Otherparameters: η1 = 8, η3 = 0.3, tend = 50. . . . . . . . . . . . . 90Figure 6.10 The power spectral densities of ψns computed using 500 re-alizations at (a) η2 = 3. (b) η2 = 2.84. (c) η2 = 2.7. Otherparameters: η1 = 8, η3 = 0.3, tend = 50. . . . . . . . . . . . . 91Figure 6.11 The coherence resonance β as a function of η2 when η1 =8,η3 = 0.3, the noise amplitude σ = 0.01. . . . . . . . . . . . 92xviiAcknowledgmentsI would like to thank all the guidance, support, and advice from my supervisorRachel Kuske. Her guidance helped me in all the time. I would also like thank mysecond reader Wayne Nagata, for giving me invaluably constructive criticism.xviiiChapter 1Introduction1.1 Bifurcation, tipping and hysteresisThe qualitative change in the behavior of ocean is often related to the bifurcation,tipping and hysteresis in the model. In brief, bifurcation refers to the qualitativechange of equilibrium as parameters vary. If the system does not follow the sameequilibrium branch as the parameters vary back and forth, it displays hysteresis.Tipping, as a result of bifurcation and hysteresis, refers to the sudden jump ofsteady state under time-dependent parameters and noise.1Figure 1.1: An example of bifurcation diagram generated by the static pa-rameter η2. η2 is the input parameter andΨ is the system response. Thesolid line is the stable steady state, and the dotted line is the unstablesteady state. See the Fig. 2.1 for detail.The Fig. 1.1 shows an example of bifurcation in a static system. In this figure,Ψ is a real-valued system response over the time t, η2 is the bifurcation parameter,and the blue line is the equilibrium as a function of η2, where the solid line is thestable solution and the dotted line is the unstable solution. For a static dynamicalsystem, a bifurcation point is where a small change in the input parameter value ofthe system induces a sudden qualitative or topological change in its behaviour[25].In our example, the upper and lower branch shows the stable equilibrium, whilethe middle dotted equilibrium is unstable. When η2 is at either L1 or L2, a littlechange in the equilibrium point structure of the underlying system will lead to asudden change in Ψ. Therefore, L1 and L2 are two bifurcations in this system.2Figure 1.2: Two examples of tipping. The left figure is generated withoutnoise and the right figure is generated with noise. η2 is the input param-eter and Ψ is the system response. The black dotted line is the equilib-rium of the underlying system (same as Fig. 1.1). The red line describesthe system when η2 is slowly increasing. See the Fig. 2.6 for detail.A dynamical system is a system where one or more functions describe the timedependence of some variables. When the parameter is time-dependent and slowlyincreasing or decreasing, the sudden change that occurs when approaching the the-oretical bifurcation point is called tipping. Tipping describes the irreversible abrupttransition in the dynamical behavior as the system moves into a different state[18].There are various types of tipping points. For example, two types are often de-scribed as B-tipping which is associated with bifurcations, and as N-tipping whichis caused by noise[2]. Here noise represents the small fluctuations in the evolutionof weather, and we usually simplify it as a random forcing in our climate models.It can be caused by forcing from land vegetation, ocean or seasonal change of cli-mate. In this thesis the tipping point is defines as a critical threshold where a smalldisturbance, either from outside noise or inside forcing, can easily alter the systemto a different state. In Fig. 1.2, Ψ is the system response and η2 is slowly increas-ing, and the red line shows how Ψ changes with η2. As η2 increases, Ψ follows thestable branch until η2 is close to L1. And because the middle part between L1 andL2 is unstable, ψ escape to the lower branch after η2 reaches L1. In the left figure,the tipping point and the bifurcation point are coincident. After adding with noise,in the right figure the tipping point is a little delayed but still close to the theoretical3bifurcation point. In many settings, tipping is associated with both bifurcations andnoise like in our right figure.Figure 1.3: One example of hysteresis behavior. Left: increasing η2. Right:decreasing η2. The red line describes how Ψ changes with η2. See theFig. 2.5 for detail.However, if we make the transition back, and find out that we can not make thesame transition at the same value of η2, then this type of behavior is called hys-teresis. Fig. 1.3 shows an example of hysteresis. The two stable steady states areseparated by some input parameter values η2 ∈ [L1,L2]. Because of this bi-stableregion, the solution switches to another stable state after tipping. Hysteresis refersto a system with two stable states that occur for the same input parameter (likeΨ at any value of η2 ∈ [L2,L1]) [27]. It is always connected with the formationof hysteresis loop: if there are two different transitions in a system that happenwhen increasing or decreasing the input power, then there is a loop between them.This interpretation emphasizes that hysteresis highly depends on its previous state,and the loop is a specific consequence when the input parameter varies back andforth[22].1.2 The Stommel modelThe meridional overturning circulation is a crucial component of the ocean circu-lation. This circulation has long been considered as a result of a variety of factorsat different latitudes in deep water across the ocean. In 1961, Stommel proposed4a model based on this (Fig. 1.4). He assumed that the overturning rate was pro-portional to density difference, and there were two critical thresholds which wouldlead to some abrupt changes inside the model. Since then, this model has beenwidely used to investigate the effect of two surface fluxes, freshwater flux and heatflux, on large-scale ocean circulation in the North Atlantic[32]. In this paper, weinvestigate hysteresis in the Atlantic circulation based on the Stommel model.Figure 1.4: two-box Stommel modelThe Stommel model describes a system with the following balance equations [32].V2dT2dt=CT2 (Ts2 −T2)+ |Ψ|(T1−T2), (1.1)V1dT1dt=CT1 (Ts1 −T1)+ |Ψ|(T2−T1), (1.2)V2dS2dt=CS2(Ss2−S2)+ |Ψ|(S1−S2), (1.3)V1dS1dt=CS1(Ss1−S1)+ |Ψ|(S2−S1). (1.4)where V1 and V2 are the volume of two boxes containing water with temperatureand salinity (T1,S1) and (T2,S2). Their subscripts 1 and 2 mean the low latitudesand high latitudes region respectively. T s and Ss are prescribed surface tempera-5ture. Due to the density gradients, heat and salt exchange in the two boxes withrelaxation coefficients CT and CS. And Ψ is the flow rate directed from high to lowpressure and is assumed linearly related to the density difference ρp−ρe betweenthe two boxes in the Stommel model:Ψ= γ · ρp−ρeρ0(1.5)andρ = ρ0(1−αT (T −T0)+αS(S−S0)) (1.6)where ρ0 is a reference density, αT and αS are the thermal expansion coefficient andsaline contraction coefficient respectively, and γ is a hydraulic constant. Consider-ing the temperature and salinity in different latitudes, we always have T1−T2 > 0and S1 − S2 > 0. If we assume the relaxation coefficients CT and CS are pro-portional to volume, then this proportion can be written as CT1V1=CT2V2≡ RT andCS1V1=CS2V2≡ RS. We set T = T1−T2, S = S1−S2, Ψ= T −S and then getdTdt= RT (T s1 −T s2 −T )−|Ψ| ·T · (1V1+1V2)dSdt= RS · (Ss1−Ss2−S)−|Ψ| ·S ·V1+V2V1 ·V2(1.7)In order to simplify, we scale time t with 1/RT , temperature T with V1V2RT/(γαT (V1+V2)), salinity S with V1V2RT/(γαS(V1+V2)), and flow rate Ψ with V1V2RT/((V1+V2)). Then Eq (1.7) turns into:dTdt=(T s1 −T s2 ) · γ ·αT (V1+V2)V1V2RT−T −|Ψ|TdSdt=(Ss1−Ss2) · γ ·αS(V1+V2)RSV1V2R2T− RSRTS−|Ψ|S(1.8)Then we define three parameters[8]:• η1—the strength of thermal forcing.– η1 = (T s1 −T s2 )γαT (V1+V2)/(V1V2RT )• η2—the strength of freshwater forcing6– η2 = (Ss1−Ss2)γαS(V1+V2)RS/(V1V2R2T )• η3—the ratio of the thermal and freshwater restoring time scales– η3 = RS/RTThe dimensionless equations now are:dTdt= η1−T (1+ |T −S|),dSdt= η2−S(η3+ |T −S|),(1.9)When we analyze ocean current, because Ψ (flow rate) is the indicator that we caneasily observe the hysteresis behavior, so we sometimes rewrite the system into Tand Ψ, where Ψ= T −SdTdt= η1−T (1+ |Ψ|)dΨdt= η1−η2−Ψ|Ψ|−T +(T −Ψ)η3(1.10)1.3 Application—tipping in the Atlantic oceanThe Atlantic Ocean circulation is one of the most critical circulation systems onearth. It carries the warmth and salinity from the low latitudes to the high latitudes.There have been several overturnings in the Atlantic Ocean circulation over thepast 2.6 million years[31]. And these overturnings result to some rapid climateshifts and have great global impacts. Therefore it is imperative to understand thedynamics of the Atlantic Ocean circulation. In this thesis, we study this problemby an analysis of the Stommel model.7Figure 1.5: The bifurcation diagram of the Stommel model when η1 = 3 andη3 = 0.3. The solid blue line is the stable solution and the dotted blackline is the unstable solution. L1 marks the point of the static non-smoothbifurcation and L2 marks the static smooth bifurcation.Fig. 1.5 shows the bifurcation diagram of the Stommel model. According tothe figure, there are two types of bifurcations in this model, one is a smooth bifur-cation and the other is a non-smooth bifurcation. After adding noise and slowlyvarying η2, tipping is induced and the main goal in this thesis is to analyze how theforcing parameter η2 and the noise lead to the transitions in the system, and thenfind ways to obtain the early warning signal before tipping. Having a better under-standing about the bifurcations in the Stommel model will provide us a perspectivefor similar types of transitions, such as the shifts in the Atlantic Ocean circulation.8Chapter 2Background on the StommelModelIn this chapter, we provide the general mathematical formulation of the dynamicsand the parameterization of the Stommel model. Considering the intricacy of atwo dimensional model, we also introduce a one dimensional hysteresis modelwith similar stability properties to the Stommel model.2.1 The static caseThe static case refers to the model with parameters that are constant, so that we cancalculate the system in equilibrium. In this section we study the static case of theStommel model, and analyze the two bifurcations in it.2.1.1 Bifurcations in the Stommel modelIn this subsection we are going to study the dynamical properties of the Stommelmodel, including the two types of bifurcations, the hysteresis behavior, and thenoise induced tipping.9Figure 2.1: The bifurcation diagram of the Stommel model when η1 = 3 andη3 = 0.3. The solid blue line is the stable solution and the dotted blackline is the unstable solution.Quantitatively, the definition of the bifurcation is when a slight variation inthe parameter value (η2 in the Stommel model) of a system causes a qualitativechange in system. Considering the process (1.10), in Fig. 2.1, the solid blue lineis its stable steady solution and the dotted black line is the unstable steady solu-tion. The bifurcation at L2 is smooth and the bifurcation at L1 is non-smooth. Forη2 < L1, only the upper branch is stable, which we call the thermally dominatedstate, whereas for η2 > L2, only the lower branch is stable which is called the salin-ity dominated state. On the branch that connects the solutions at L1 and L2 (dottedline in Fig. 2.1), the steady state is unstable. The system when L1 < η2 < L2 mayevolve to either the upper or lower branch depending on its previous state.102.1.2 Location of the smooth saddle node bifurcationTo categorize the smooth bifurcation at L2, we first consider the steady state on theupper branch. Writing T as a function of Ψ, and Ψ as a function of T,η2, we have:Ψ˙= f (Ψ,η2) = η1−η2−Ψ2−T (Ψ)+η3(T (Ψ)−Ψ)T˙ = g(T,η2) = η1−T (1+Ψ)If one assumes the smooth bifurcation is at (Ψ∗, η∗2 ), it can be inferred that f (Ψ∗,η∗2 )=0, g(T ∗,η∗2 ) = 0. Then we do Taylor expansions for f (·) and g(·) at (Ψ∗, η∗2 ) and(T ∗, η∗2 ) respectively, we get:Ψ˙=−(η2−η∗2 )+(−1+η3−12TΨΨ(Ψ∗))(Ψ−Ψ∗)2T˙ = (−1−Ψ∗−T ∗ΨT (T ∗))(T −T ∗)+ 12(−2ΨT (T∗)−T ∗ΨT T (T ∗))(T −T ∗)2Next, let β =η2−η∗2 and P=Ψ−Ψ∗, scaling P and β with A=(−1+ η3−12 TΨΨ(Ψ∗)),we haveP˙ =−Aβ +P2When T˙ = 0, we haveT =η11+ΨSo thatTΨΨ(Ψ∗) =η11+Ψ∗The smooth bifurcation Ψ∗ is positive as shown in Fig. 2.1, and we have η3 < 1(η3 is the ratio of the thermal and freshwater restoring time scales, or refer to thevalue in the Table (2.1)), so that A 6= 0. Now we can conclude that L2 is a saddlenode bifurcation based on the form of P˙.2.1.3 Behavior near the equilibriumTo derive the dynamics near the steady solution shown in Fig. 2.1, we need tolinearize the following equations:11 T˙ = η1−T (1+ |T −S|)Ψ˙= η1−η2−Ψ|Ψ|−T +(T −Ψ)η3 (2.1)For Ψ> 0 T˙ = η1−T (1+Ψ)Ψ˙= η1−η2−Ψ2−T +(T −Ψ)η3 (2.2)If we suppose that the steady solutions are Teq and Ψeq satisfying:T˙ = 0→ Teq = η11+ΨeqΨ˙= 0→ η1−η2−Ψ2eq−η11+Ψeq+η3(η11+Ψeq−Ψeq) = 0 Ψeq > 0and define u = T − Teq, ψ = Ψ−Ψeq, then we may linearize the system aboutequilibrium: [u˙ψ˙]=[−1−Ψeq −Teq−1+η3 −2Ψeq−η3][uψ]The eigenvalues of the above matrix are:λ =12(−η3−3Ψeq−1±√η23 −4η3Teq+2η3Ψeq−2η3+4Teq+Ψ2eq−2Ψeq+1)To simplify, we would like to set a reasonable value to η3. Recall that η3 is the ratioof the thermal and freshwater restoring time scales, so that η3 is positive. Besides,thermal restoring time is longer than freshwater[35], thus η3 ∈ (0,1). In Fig. 2.2we test four different values of η3. As η3 increases toward 1, the smooth saddlenode bifurcation disappears gradually.12Figure 2.2: The bifurcation diagram of the Stommel model when η1 = 3 andη3 = 0.1,0.3,0.5,0.8 respectively.Here we use the same values of parameters as in the chapter 6 in [8]. Forη3 = 0.3 and Ψ> 0, after simplificationλ =12(−η3−3Ψeq−1±√0.49+8.41+Veq+V 2eq−1.4Veq)The system always has two real eigenvalues. Also, forΨ ∈[0, 0.39] (approximate),one of the eigenvalues is positive, indicating that the solution is not stable as wehave shown in Fig. 2.1.Figure 2.3: The steady state of Eq (2.2). The dotted black state is unstable,and the red state is stable. η1 = 3, η3 = 0.3.13For Ψ< 0 T˙ = η1−T (1−Ψ)Ψ˙= η1−η2+Ψ2−T +(T −Ψ)η3 (2.3)On this lower branch, Teq and Ψeq satisfyT˙ = 0→ Teq = η11−VeqΨ˙= 0→ η1−η2+Ψ2eq−η11−Ψeq +η3(η11−Ψeq −Ψeq) = 0 Ψeq < 0After linearizing the system (2.3) about the equilibrium Teq and Ψeq, we get[u˙ψ˙]=[−1+Ψeq Teq−1+η3 2Ψeq−η3][uψ](2.4)The eigenvalues of the above matrix are:λ =12(−η3+3Ψeq−1)±12√η23 +4η3Teq−2η3Ψeq−2η3−4Teq+Ψ2eq+2Ψeq+1) (2.5)When η3 = 0.3 and Ψ < 0, the steady solution of the system is always stable.Since η2 represents the strength of the freshwater forcing, we only take positivevalues of η2 into account. When η2 ∈ [0.9,8.32](approximation), the matrix (2.4)has two complex eigenvalues λ with negative real parts, and when η2 ∈ [8.32,+∞],the matrix (2.4) has two negative real eigenvalues λ . The red solid line in Fig. 2.4is the stable steady state.14Figure 2.4: The steady state of Eq (2.3). The red solid line is the stable solu-tion. η1 = 3, η3 = 0.3.2.2 Varying parameterDue to the feedback between the flow and salinity transportation, η2 changes withtime, which results in the change of system state. Fig. 2.5 shows how the systemalters its state under different varying direction of η2. When η2 increases fromthe smaller value, it follows and finishes the upper branch, jumps to the lowerbranch through the smooth bifurcation. While if η2 decreases from the larger value,it follows the lower branch until the non-smooth bifurcation point. These twotransitions for increasing η2 and decreasing η2 model the origin of some large andabrupt fluctuations in the high-latitude climate in North Atlantic.15Figure 2.5: The bifurcation diagram when η2 is slowing varying. Left: η2increases. Right: η2 decreases. η1 = 3, η3 = 0.3.2.3 Possible noise sourcesHowever, it is not always the case that the system we study follows the dynamicsin Fig. 2.5. If we review how the coefficients are defined in section (1.2):• η1 = (T s1 −T s2 )γαT (V1+V2)/(V1V2RT )• η2 = (Ss1−Ss2)γαS(V1+V2)RS/(V1V2R2T )the variables that constitute them in fact change with the ocean conditions. Ta-ble (2.1) describes those parameters based on some references [21][35]. The hy-draulic constant γ shows the ease with which a fluid (usually water) can move, andin this case, γ is chosen to obtain a reasonable value of the THC strength. If weconsider the seasonal changes that are related to the changes in the strength andstability properties of the currents, we may conclude that γ also has the seasonalfluctuations[33]. Therefore we put noise in η1 and η2, and since they are additive,we use σWt to represent them. As a result, we have the Stommel model Eq (2.6)with additive noise:dTdt= η1−T (1+ |T −S|)+σ1Wt ,dSdt= η2−S(η3+ |T −S|)+σ2Wt .(2.6)16Table 2.1: Parameters of the box modelParameter Symbol ValueHydraulic Constant γ 1×109m3s−1Thermal Expansion Coefficients αT 1.5×10−4K−1Saline Relaxation Coefficient CS 1.9×106m3s−1Volume V1 1017m3Volume V2 1017m3C over V , is a constant RT 3.17×10−11s−1C over S, is a constant RS 1.9×10−11s−1Prescribed temperature for polar region T s2 0◦CPrescribed temperature for equatorial region T s1 30◦CPrescribed salinity for equatorial region Ss1 35.6psuPrescribed salinity for polar region Ss2 34.7psuThe fluctuations may cause noised-induced tipping [3] and add uncertainties toour prediction. To better illustrate, one may check the time series in Fig. 2.6 thatis generated with noise and a slowly increasing η2. In this case, tipping occursbeyond the static bifurcation value, but it can occur earlier for other realizations ofthe noise. Therefore, how to extract useful information from the noisy data, so thatwe can predict exact tipping time is the focus in this thesis.17Figure 2.6: One realization generated from Eq (2.6). η2 is slowing increas-ing. η1 = 3, η3 = 0.3, dt = 0.1, ε = 0.01, σ1 = σ2 = 0.1.Thompson and Sieber[30] considered a model with additive noise of this type:dx = f (x,λ (t))dt+g(x)dW (2.7)where W (t) represents Brownian motion.2.4 Simplified 1D Hysteresis modelIn order to study the transitions in Eq (2.6), we introduce a 1D hysteresis model (2.8)with a smooth and a non-smooth bifurcation. Even though the form is simple, this1D model has similar properties compared with the 2D Stommel model (1.10). x˙ =−µ+2|x|− x|x|µ˙ = ε (2.8)Fig. 2.7 shows the bifurcation diagram of the static 1D hysteresis model (2.8).It has the same types of bifurcations as the Stommel 2D model.18Figure 2.7: The bifurcation diagram of the 1-D hysteresis model (2.8). Thesolid line is the stable steady state and the dotted line is the unstablesteady state.2.4.1 Bifurcations in the 1D modelThe smooth fold bifurcationThe hysteresis model (2.8) is composed of two processes. The one that is usedwhen x > 0 is x˙ = (−µ + 1)− (x− 1)2, and the one that is used when x < 0 isx˙ = (−µ−1)+(x−1)2. Fig. 2.8 shows the steady state of the process x˙ = (−µ+1)− (x−1)2 used when x > 0. The red line is the stable solution and the black lineis the unstable solution. The dotted line is when x < 0 which should be replacedby another process x˙ = (−µ−1)+(x−1)2.Figure 2.8: The steady state of process x˙ = (−µ+1)− (x−1)2.19If we set µ1 = −µ + 1, X1 = x− 1, then the process on the upper branch isequivalent toX˙1 = µ1−X21 (2.9)which is a saddle-node form.The non-smooth fold bifurcationWhen x< 0, the 1D hysteresis model turns into x˙=(−µ−1)+(x−1)2. The steadystates are shown in Fig. 2.9. The dotted line does not show up in the 1D modelbecause this figure is only applicable for x < 0. The dotted red line is the unstablesteady solution and the solid red line is the stable steady solution. The solid red linein Fig. 2.9 and the solid lines in Fig. 2.8 together make up the bifurcation diagramfor Eq (2.8).Figure 2.9: The steady state of process x˙ = (−µ−1)+(x−1)220Chapter 3Noise-induced tipping near asmooth saddle node bifurcationThe smooth bifurcation in the Stommel model is proved to be a saddle node bi-furcation in Section (2.1.2). Fig. 3.1 shows the bifurcation diagram of the saddlenode normal form x˙ = µ+ x2. Once noise is added to this form, tipping will be in-duced. We study the properties of this system just before tipping and hence providea possible method to predict tipping in the Stommel model (2.6).Figure 3.1: Bifurcation diagram of the saddle node form x˙ = µ+ x2.In 2012, Sieber introduced and tested several approaches, including linear de-21cay rate, potential well and skewness, to extract tipping information from the twoprocesses (3.1) and (3.2) below [30].dx = [µ− x2]dt+σdWtdµ =−εdt(3.1)dx =−2√µxdt+σdWtdµ =−εdt(3.2)Here x is the response, µ is a parameter with drifting rate −ε , t is time, Wt isthe Brownian motion with zero mean and variance t so that σ controls the noiseamplitude. Process (3.1) is a saddle-node normal form with a drifting parameter µand additive noise. Process (3.2) is the linearized form of (3.1). They both have anidentical linear decay rate 2√µ , but are different in nonlinear terms: process (3.1)has quadratic term x2 but process (3.2) does not.3.1 Time series analysis methodsIn this section, we introduce and provide the mathematical explanation of the threeapproaches — linear decay rate, potential well and skewness. Then these meth-ods are applied in the processes (3.1) and (3.2) to derive the linear and nonlinearproperties in them.3.1.1 Gaussian kernel detrendingGaussian kernel detrending is the first step to process our data. It helps remove thetrend of the data, so that we can focus on the behavior of the fluctuations.22Figure 3.2: In both plots, the black curve represents the time series, the redline is the theoretical equilibrium, and the blue curve is the result of fil-tering the time series with a Gaussian kernel. The left plot is obtainedfrom the process (3.1), and the right plot is from the process (3.2). Pa-rameters: µ(0) = 2, ∆t = 0.05, σ = 0.35, ε = 0.009, bandwidth d = 8.In Fig. 3.2, the black curve represents our data generated within model (3.2)and model (3.1), the red line is the theoretical equilibrium of each model: the de-terministic saddle-node form x˙ = µ − x2 has the stable equilibrium xeq =√u andthe unstable equilibrium xeq = −√µ , the linear time series x˙ = −2√µx has thestable equilibrium xeq = 0. We may notice the equilibrium value (red) of x frommodel (3.1) changes over time. Thus if we use the original time series (black) inthe Fig. 3.2 directly, the varying mean might cause some kind of distortion. Henceit is necessary to remove this trend so that one can focus on what remains in the lefttime series in Fig. 3.2. When µ increases towards the critical value µc (eg. µc = 0in nonlinear process (3.1)), the noise is likely to cause tipping [30] : as we see inthe left Fig. 3.2, x suddenly drops to a much smaller value when µ approaches 0.If one can detect an increased probability for tipping from the detrended data, thenit is possible to predict tipping.There are various ways to remove the trend. Lenton[17] proposed fitting theoriginal data with a scheme based on a Gaussian kernel. A Gaussian kernel is akernel with shape of a normal distribution. It is widely used and makes a gooddefault kernel. Regardless of the Gaussian width, one property is that Gaussiankernel is universal. It can be used to approximate a bounded continuous function23up to arbitrary precision needed (in the sense of ‖ · ‖∞). We use a Gaussian kernelas followsGk(t) =1√2pidexp(−12(t− k∆t)2d2) (3.3)N is the total data points, k is the index of each point, k = 1,2, ...,N, and d is thebandwidth which controls the smoothness or roughness of a density estimate. Sinceµ is varying on a slow time scale compared to that of the time series (µ could bethe rate of change over thousands of years of an ocean, while the data is collectedseveral times per year), then to capture the moving trend on a slow time scale andalso reflect the variation of the data, the bandwidth needs to satisfykdri f t 1d kcrit . (3.4)where kdri f t is the averaged drifting rate of a parameter, which is εt in this process.The rate kcrit is the the rate with which a small disturbance in state space returnsto the steady state, or to say, the linear decay rate. Then to remove the trend, wesubtract the kernel-weighted average X(k∆t) (blue line in Fig. 3.2) from originaldata and get the detrended data. Each X(k∆t) is estimated using all the data: thecloser a data point to the kth point, the more weight that data point has on X(k∆t).In this way we get a smooth moving average,X(k∆t)k=1,2,...,N =∑Ni=1 Gk(i∆t)xi∑Ni=1 Gk(i∆t)(3.5)The detrended time series x˜ isx˜k = xk−X(k∆t),k = 1,2, ...,N (3.6)x˜k fluctuates around zero on the time scale of u according to the Eq (3.4): it can beconsidered as stationary on time scales faster than εt, which is 1kdri f t [14].3.1.2 Estimate of the linear decay rateIn Chapter 2, we have discussed that when µ increases toward µc, the noise islikely to cause tipping. This behavior can be characterized through the linear de-24cay rate which captures the strength of attraction for the attracting branch. Thedecay rate depicts how fast the state relaxes back to the stable branch after a smalldisturbance. The smaller the decay rate, the longer the system takes to return to theattracting branch. For µ far from the bifurcation, x can return back to the equilib-rium quickly. For µ close to the bifurcation, the system takes longer time to return,or even jumps away from the current state. In our example (3.1), x is more likely toescape when µ is close to µc. Using the detrended variable x˜ = x−X derived fromEq (3.6), the next step is to calculate the decay rate toward the trend X and its timedependence. Here we introduce two schemes to compute the decay rate.The first scheme is based on the autocorrelation function (ACF) and we call thelinear decay rate kACF . It computes the correlation between successive elementswithin a window sliding along the time series[5]. The dynamics of the 1D timeseries can be interpreted as a one-step autocorrelation function (ACF):x˜k+1 = α x˜k +σwk, (3.7)where σ controls the additive noise strength added to the system and wk is Gaussiannoise. The linear decay rate kACF of the system is related to α viaα = exp(−kACF∆t). (3.8)Hence when µ is approaching the bifurcation point µc, kACF vanishes and αgoes to one. This was introduced by Held and Kleinen[15] when analyzing climatetime series. Because the drifting rate of our parameter is slow (ε 1), α can beseen as a constant on a time scale slower than εt, and that can be true if we usea window size w = 2m+ 1 satisfying w < 1∆tkdri f t when simulating α [30]. Forexample, we can determine αk by fitting points inside the kth window to Eq (3.7),and αk is the coefficient in the linear relationship: x˜ j+1 = αkx˜ j, j = k−m, ...,k+m− 1. The Matlab package regress is used to do multiple linear regression inthis thesis: for a given data set, it determines the best fit via the least squares.This method supposes the form of function as y = f (x,α), and in our problemf (x,α) = αx. For each point xi+1, the error is ri = xi+1− f (xi,α). The goal is to25find a value of α that minimizes the sum of squares of errors ∑r2i : that is when∂ ∑r2i∂α= 0,ri = xi+1− f (xi,α) (3.9)and in our problem, we estimate αk in the kth window via minimizingk+m−1∑i=k−m(xi+1−αxi)2In numerical calculation, it should be treated with caution when choosing awindow length. Fig. 3.3 displays the linear decay rate of four cases computedusing four window lengths. Window length affects two key aspects - temporalresolution and coefficient quality. A proper window size should lead to a generaltrend in the coefficient without large fluctuations. When the window size is small,we have more windows, and each window reflects greater temporal resolution be-cause the time period over these points is short. However, the result also has largerfluctuations and lower coefficient quality, which refers to the smoothness of thecurve. Large fluctuations make it hard to tell the real value of the linear decayrate since they are mainly noise. As we see for the cases of small window lengthN/4 and N/5 in Fig. 3.3, we have high resolution but the fluctuation on the leftside is too large. On the other hand, for large window size, we then have fewerwindows across the data, and thus for each window, there is a long time interval.To calculate the linear decay rate at a specific time, we use data that is far beforethat time or much later. As a result, the linear decay rate can not reflect the latestchange as it is diluted over a long time period. Hence the temporal quality is poorif the window size is too large, although the coefficient is smoother. As shown inFig. 3.3, when using long window length N/2 to compute the linear decay rate, themethod uses half the data and gives less temporal resolution. The black and the redlines are smooth but not able to capture the information locally. A suitable windowlength should not be too large so that it can capture the recent variation, and alsonot too narrow to avoid noise. Therefore, in this section we choose a window sizeas N/3. It is also noticeable that after applying a window w, the output length oflinear decay rate is [w/2,N−w/2]26Figure 3.3: The linear decay rate for Eq (3.1) with different window lengths.The black line is the kACF , the red line represents kU as defined below.We try four different window lengths, and N is the total number of ourdata. Top left: the decay rate with window length N/2, top right: N/3,bottom left: N/4, bottom right: N/5. Other parameters: µ(0) = 2,∆t = 0.05, σ = 0.35, ε = 0.009, bandwidth d = 8.Besides ACF, one can also estimate the decay rate via approximation for theeffective potential well and we call this decay rate kU . The potential well U is thethe energy landscape surrounding a local minimum of potential energy. The poten-tial well and system response can be related by a stochastic differential equationx˙(t)=−U ′(x)+σdWt . For example, in process (3.2), U =√ux2. For process (3.1),where x˙ = µ− x2+σdWt , U = 13 x3−µx.For the detrended time series x˜, the minimum of the potential well gives thelocation of the equilibrium or the steady state: x˜ = 0, and we have ∂xU(0) = 0. Ifone assumes the deterministic dynamics of x˜ is essentially that of an overdampedparticle moving in a potential well U(x˜), then the stochastic expression is:dx˜ =−∂xU(x˜)dt+σdWt (3.10)We linearize ∂xU at the equilibrium x˜ = 0 and since ∂xU(0) is equal to 0 after27detrending the time series x, we have ˙˜x = −x˜ ·Uxx(0)+σdWt . Hence the lineardecay rate kU equals the second derivative ∂xxU(0). Moreover, the shape of U isrelated to the stationary density. According to the Fokker-Planck equation, for astochastic differential equationdXt = ω(Xt , t)dt+σ(Xt , t)dWtwith drift ω(Xt , t) and diffusion coefficient D(Xt , t) = σ2(Xt , t)/2, the Fokker-Planck equation for the probability density p(x, t) of the random variable Xt is∂∂ tp(x, t) =− ∂∂x[ω(x, t)p(x, t)]+∂ 2∂x2[D(x, t)p(x, t)] (3.11)The stationary solution is given when ∂ p∂ t = 0. It can be inferred that for pro-cess (3.10) we have12∂xx p(x˜) = ∂x[−σ−2∂xU(x˜)p(x˜)], (3.12)which gives, integrated once,12∂x p(x˜) =−σ−2∂xU(x˜)p(x˜)+σ−2c, (3.13)where c is the constant of integration. Since the probability density becomes zeroat ±∞, c = limx˜→±∞ ∂xU(x˜)p(x˜). Here we expand −σ−2∂xU again at the bottomof the well x˜ = 0, remembering that ∂xU(0) is equal to 0 and the decay rate kU isequal to ∂xxU(0):12∂x˜ p(x˜) =−σ−2x˜ ·∂xxU(0)p(x˜)+σ−2c=−σ−2x˜ · kU · p(x˜)+σ−2c=−k′U x˜ · p(x˜)+σ−2cwhere kU = k′Uσ2. In order to determine the coefficient k′U in the equation, wefirst need to obtain the empirical density pemp(x˜) in every sliding window alongthe time series. That can be estimated using the Matlab package ksdensity. ksden-28sity is based on a normal kernel function[26]. The kernel density estimator is theestimated probability density function of a random variable. For any real values ofx˜, the kernel density estimator’s formula is given bypemp(x˜) =1NhN∑i=1K(x˜− x˜ih) (3.14)where N is the number of points from the unknown distribution, h is the defaultbandwidth that is the optimal for normal densities[4].h = σˆ(43N)15Here σˆ is the standard deviation of the sample data. K(·) is the kernel smoothingfunction. In this thesis, we use the standard normal distribution :K(x) =1√2pie−12 x2Then to get k′U , we use relation (3.15) to find k′U from pemp12∂x pemp(x˜) =−k′U x˜pemp(x˜)+ cemp (3.15)via nonlinear curve-fitting using the Matlab package lsqcurvefit [7], where k′U =σ−2kU and cemp = σ−2c. In each sliding window, we use pemp(x˜) estimated fromEq (3.15) to obtain the k′U . For function y = F(c,x), where in this context y isgiven by 12∂x pemp(x˜) and F = −k′U x˜pemp(x˜)+ cemp, lsqcurvefit finds coefficient cby solving the problem in the kth window (length 2m+1)minc‖ F(k′U ,x)−12∂x pemp ‖22= minck+m−1∑i=k−m(F(k′U ,xi)−12∂x pemp,i)2 (3.16)via the trust-region-reflective least squares algorithm[6]. The sketch of trust-region-reflective least squares is:First, we formulate a sub-problem min{12 sT Hs+sT g such that ‖ s ‖≤∆}, whereg is the gradient of (F − y) at current c value, H is the Hessian matrix that is thesymmetric matrix of the second derivatives of (F−y), ∆ is a positive scalar that de-29scribes the trust-region dimension, and ‖ · ‖ is the 2-norm. Our goal is to solve fora trial step s in the neighborhood that has ∑i(F(c+ s,xi)− yi)2 < ∑i(F(c,xi)− yi)2.More´ and Jorge[24] introduced a algorithm that involves the computation of a fulleigensystem and a Newton process to the secular equation1∆− 1‖ s ‖ = 0Then if ∑i(F(c+ s,xi)− yi)2 < ∑i(F(c,xi)− yi)2, update c = c+ s, and adjust ∆according to standard rule: if the current iteration makes a satisfactory reduction,we may exploit our model more in the next iteration by setting a larger ∆. If weonly achieve a limited improvement after the current iteration, the radius of thetrust-region then should not have any increase. If the trial step s is not accepted, wedecrease ∆ and try again. We keep repeating the above steps until convergence andtherefore get the coefficients that minimize the objective function. To sum up, ineach sliding window, we estimate the empirical density pemp(x˜) by the kernel den-sity estimator Eq (3.14). And then the linear decay rate k′U for that sliding windowis obtained by fitting it to Eq (3.15).In addition to using Eq (3.15) directly, one can also use the solution form ofEq (3.15)pemp = 2cempe−k′U x2∫ek′U x2dx+√k′Upie−k′U x2(3.17)in lsqcurvefit to fit k′U to pemp. Then F is the right side of Eq (3.17), c is kU andcemp, and y is pemp. In this paper, what we choose to use is Eq (3.15). We derivek′U through curve fitting, and the exact linear decay rate kU can then be determinedif we have the noise amplitude σ . kU and kACF measure the same thing and shouldequal to each other. In the real world, it is usually hard to obtain the noise ampli-tude σ . If kU = kACF , once we have kACF and k′U , the noise amplitude σ can beestimated via σ2emp = kACF/kU .The red line in Fig. 3.4 depicts how α changes with time and the black lineshows the estimation of kU from the empirical density. The black and the red linesare quantitatively close to the theoretical line, and are strongly correlated with each30other. When µ is close to the bifurcation point, the system is more likely to tip withadditive noise. Thus α in Eq (3.8) approaches 1 and the linear decay rate decreasesto 0.However, as to the question about the timing of tipping, the result of the linearanalysis of decay rate is not straightforward. The two time series have identicalproperties as seen in Fig. 3.4, while the tipping is widely different in these twocases. For the linear process, its tipping occurs much later than the nonlinear pro-cess, but this could not be deduced from Fig. 3.4. The trend of kACF and kU is anindicator for increased likelihood of tipping while it neglects the role played by thenonlinear term.Figure 3.4: The right one is the result of the linear decay rate coming fromthe process (3.1) and the left one is coming from the process (3.2). Theblack curve represents kU and the red curve is kACF , the smooth blueline is the theoretical decay rate exp(−2√µ∆t) calculated from the lin-earized deterministic equation. The window length is w = N/3.3.1.3 Nonlinearity—information from the potential wellWe seek supplementary information about the nonlinear features. In the climatesystem, tipping is a feature of nonlinear dynamics. Analysis in absence of theleading nonlinear terms can not provide accurate prediction for early escape [30].We continue our analysis of the potential well, which was used to obtain kU .We already know from Section (3.1.2) that a stochastic differential equationcan be written in the formx˙(t) =−U ′(x)+σdWt31Hence process (3.1) can be rewritten as:x˙(t) =−Ux+σdWt ,whereUx = x2−µ(t) (3.18)U =13x3−µ(t)xAt each discrete time, the potential well U(x) has a shape as quadratic.The process (3.2) can be rewritten as:x˙(t) =−Ux+σdWtwhere Ux = 2√µ(t)x (3.19)U =√µ(t)x2At each discrete time, the potential well of the linear process has a shape as parabola.Raymond[29] states that the stable and unstable equilibrium can be marked bypotential curves that are concave or convex for a general potential function U(x).Furthermore, in the immediate neighborhood of a stable equilibrium, say it is atx = a, U(x) can be fit quite well by a parabola:U(x) =12K(x−a)2where K = d2Udx2 |a.The form of the potential well U can also be obtained through a similar cal-culation as shown in section (3.1.2) via Fokker-Planck equation: for a stochasticdifferential equationdXt = ω(Xt , t)dt+σ(Xt , t)dWtwith drift ω(Xt , t) and diffusion coefficient D(Xt , t) = σ2(Xt , t)/2, the Fokker-32Planck equation for the probability density p(x, t) of the random variable Xt is∂∂ tp(x, t) =− ∂∂x[ω(x, t)p(x, t)]+∂ 2∂x2[D(x, t)p(x, t)]The stationary solution ∂ p∂ t = 0 is given by∂∂x[ω(x, t)p(x, t)] =∂ 2∂x2[D(x, t)p(x, t)]And in our problem, it is12∂xx p(x) = ∂x[−σ−2∂xU(x)p(x)] (3.20)Gardiner[10] gives a stationary solution of the Fokker-Planck equationp(x) = exp[−2U(x)/σ2]The empirical potential well can be estimated from the time series viaUemp =−σ2emp2log pemp. (3.21)As we see in Equation (3.15), k′U is σ−2kU where kU and σ are linear decay rateand input noise amplitude respectively, so we may replace σ2 with kACF/kU ob-tained from the autocorrelation function Eq (3.15) and the potential well Eq (3.13)respectively:σ2emp = kACF/kU (3.22)In a real world problem, it is hard to know the amplitude of noise, thereforestep (3.22) is essential. Then we extract the empirical density pemp(x) in everysliding window via the probability density estimate along the time series. That canbe done with the Matlab package ksdensity. In the kth window with length 2m+1,we take points [xi]k+mi=k−m and solve the probability density at time tk (see Eq (3.14)for details)pemp,k(x) =1Nhk+m∑i=k−mK(x− xih) (3.23)33Fig. 3.5 shows the potential well of process (3.1) and process (3.2). The timeseries x is displayed as a blue curve at the bottom of the Uemp valley. In every timewindow, we have one empirical noise amplitude as a function of x derived fromthe linear decay rate and the empirical probability density derived from Eq (3.23).Then we can compute the empirical potential well centered at that time windowusing Eq (3.21) which is shown as the grey surface. We approximate both po-tential wells by parabola. Then we study the deviation of the potential well fromits approximated parabola, and compare between the linear and nonlinear process.Here we derive the form of parabola via the Vandermonde matrix (or by the Matlabpackage polyfit):• Suppose we have M query points (a1,a2, ...,aM), N fitted values (b1,b2, ...,bM)at query points and a predetermined polynomial order n.• Form the Vandermonde matrix Van1 an−11 ... 1an2 an−12 ... 1... ... ... ...anN an−1N ... 1• Solve for coefficients P with p =V\ban1 an−11 ... 1an2 an−12 ... 1... ... ... ...anM an−1M ... 1P1P2...Pn+1=b1b2...bM (3.24)We focus on the potential well computed from the last window because it containsthe most recent data before tipping. Assuming it is the kth window, first we estimatethe probability density (Eq (3.23)) using points xk−m,xk−m+1, ...,xk+m. The outputpemp is a function of x centered at tk and is evaluated at equally-spaced points ai thatcover the range of the data. We name these points as (a1, p1),(a2, p2), ...,(aM, pM).Then the empirical potential well centered at tk can be determined by Eq (3.21).The potential well is also a function of x: (a1,U1), (a2,U2),...,(aM,UM). We sup-34pose the form of the potential well as P1a2 +P2a+P3, then P can be worked outby a21 a11 1a22 a12 1... ... ...a2M a1M 1P1P2P3=U1U2...UMThe red line in Fig. 3.5 is the parabola approximated to Uemp from the last window.The differences between the approximated parabola and the empirical potentialwell is measured by the normalized root-mean-square deviation NRMSD wheresmall values indicate small residual variance:NRMSD =√∑Mi=1(Ui− (P1a21+P2a22+P3))2√M(Umax−Umin)(3.25)The NRMSD of the linear process is 0.049 and of the nonlinear process is 0.069.The right one in Fig. 3.5 is the potential well of the linear time series. FromEq (3.19), its potential well is already a parabola, so we see it overlaps the ap-proximated parabola shown in red color. However, in the left Fig. 3.5, the potentialwell does not fit the parabola as well as in the linear process. Compared with thered parabola, the potential well is slightly above it for the larger values of x whilebeneath it for the smaller values of x. Eq (3.21) shows that the depth of the poten-tial well is proportional to the probability density. Thus the part of the empiricalpotential well that falls below parabola indicates that there are more points in thatdirection, and a tipping to this direction is pending.35Figure 3.5: The left figure is the potential well of the time series (3.1). Theright one is the linear time series (3.2). The grey surface is the potentialwell estimated by Eq (3.21), the red line is a parabola approximated tothe potential well via Eq (3.24), and the blue curve is the time series.To better illustrate the approximation to a parabola for the nonlinear process (3.1),Fig. 3.6 shows the theoretical (grey) and empirical (black) potential well and thebest fit parabola (red) for that process. When u is far from the critical point µ = 0,these three lines fit pretty well with each other and the potential well can be fitquite well by a parabola, which is consistent with the statement in Raymond[29].However, the final theoretical grey line is more shallow towards the negative valuethan the other two before it, which indicates the system is more likely to jumpto the smaller value by that time. Furthermore, because the left side of potentialwell is more shallow, when we approximate the empirical potential well with theparabola, the parabola tends to rise above the empirical and theoretical potentialwell at the negative value side and falling beneath them at the positive value side.Then we may conclude that a jump to the smaller x may be pending.36Figure 3.6: For nonlinear process (3.1): the black dotted line is the empiricalpotential well, the red solid line is the approximated parabola to the em-pirical potential well, and the grey solid line is the theoretical potentialenergy.3.1.4 SkewnessThe other method that has been widely proposed to detect nonlinearity is skew-ness. Guttal and Jayaprakash studied how bifurcation influenced skewness [11, 12].Skewness measures the asymmetry of the empirical quasi-stationary density of avariable about its mean. For example, for a Gaussian distribution, its density isevenly distributed around mean value and its skewness is 0. By contract, for anonlinear saddle node process with additive noise, the system spends more timereturning back to equilibrium when µ is near µc, which results in a highly skewedempirical probability density in the direction of bifurcation. The definition of skew-ness γ is (Matlab package skewness) :γ =∫(x−mean(x))3 p(x)dxσ3. (3.26)For our discrete time series, this is:γ =1n ∑ni=1(xi− x¯)(√1n ∑ni=1(xi− x¯)2)3By plotting the skewness as a function of time for slowly varying µ , the non-37linear behavior can be detected. Fig. 3.7 shows the estimate of skewness generatedby 100 realizations of Equation (3.1) and Equation (3.2). Each time series is rununtil N = 2000, or x =−1 (indicating tipping). We save 100 trajectories that havenot escaped by N = 2000. And then for each discrete value of µ , we use 100 pointscoming from 100 different time series to calculate the skewness. It is noticeablethat the skewness is negative in the nonlinear case and decreases when approachingthe bifurcation point. But in the linear system, γ remains close to zero, which isconsistent with the observation that in this time period, the system is less likely totip. Therefore, an increase in the magnitude of the skewness is another possibleearly warning signal.Figure 3.7: The upper plot is the skewness γ for the process (3.1) and thelower plot is for the process (3.2). Skewness is calculated from 100realizations.38Chapter 41-D model with smooth andnon-smooth bifurcationsTo extract tipping information from the Stommel model, we first apply all methodsthat have been investigated in Chapter 3 to the simplified 1D model. The sim-plified 1D model is introduced in Section 2.4. This model has the same types ofbifurcations, i.e., one smooth and one non-smooth bifurcations compared with the2D Stommel model:y˙ =−υ+2|y|− y|y|,Figure 4.1: The bifurcation diagram of the 1-D model. The dotted line is theunstable steady state. The solid line is the stable steady state.39where υ is the drifting parameter and y is the system response. Fig. 4.1 is the bifur-cation diagram of the 1-D model. The black line is the steady state of y. For υ < 0,only the upper branch is stable, whereas for υ > 1, only the lower branch is stable.On the dotted branch 0< υ < 1, the steady state is unstable, the system may evolveto either the upper branch or the lower branch. In this model the smooth bifurcationis at υ = 1, and the non-smooth bifurcation is at υ = 0.When υ is slowly varying with rate υ˙ = δ , the system alters its state differentlyfor increasing or decreasing υ . We observe hysteresis because of the differenttransition values of the bifurcations: when υ increases from the smaller value, itfollows the upper branch and then passes through the static smooth bifurcation,while if υ decreases from a larger value, it follows the lower branch and passesthrough the static non-smooth bifurcation.Figure 4.2: The hysteresis behavior in the 1-D model. The bifurcation dia-gram when υ is slowly varying. Left: when υ is increasing. Right:when υ is decreasing.The 1-D process with additive noise is y˙ =−υ+2|y|− y|y|+σWtυ˙ = δ (4.1)where σ controls the noise amplitude added to the system and Wt is Brownianmotion with mean zero and variance t. To make the 1D model comparable with themodel (3.1), we set σ = 0.35 (see Appendix (A.1) for detail).404.1 Location of the smooth saddle-node bifurcationIn this section, we apply the three methods — linear decay rate, potential well andskewness in the process (4.1) to derive the linear and nonlinear properties as υapproaches the smooth saddle-node bifurcation.4.1.1 Time series and Gaussian kernel detrendingWe first generate a time series using the process (4.1), and then extract the partbefore tipping to remove the trend.Figure 4.3: The time series coming from the process (4.1) where υ is slowlyincreasing. The red line is the equilibrium, the green line is the resultof filtering the time series with a Gaussian kernel, and the black line isthe system response ys. Parameters: σ = 0.35, δ = 0.009, ∆t = 0.05,bandwidth d = 8.We use ys to represent the system when υ varying through the smooth staticbifurcation. The black curve in the Fig. 4.3 is the time series generated by theprocess (4.1), the red line is the theoretical equilibrium and the green line is theGaussian moving average. In this realization tipping is at t = 0.95, where ys sud-denly drops dramatically. As shown in Fig. 4.3, the whole time series is fromυ =−1 until υ = 1.5, but what we should focus on is the part before tipping. Wepick out that part υ ∈ [−1,0.9], supposing the number of points inside that regionis N, then we average ys over a Gaussian kernel Gk(t)k=1,2,...,N that is introduced in41Eq (3.3), and get the smooth fitting line Ys(k∆t)k=1,2,...,N byYs(k∆t)k=1,2,..,N =∑Ni=1 Gk(i∆t)ys,i∑Ni=1 Gk(i∆t)(4.2)The moving average Ys is indicated by the green line in Fig. 4.3. After determiningthe trend Ys, the detrended time series is estimated by y˜s = ys−Ys (Fig (4.4)).Figure 4.4: The detrended data y˜s of the time series (black) in Fig. 4.3, ob-tained by subtracting the moving average y˜s (green line in Fig 4.3) fromthe original time series ys (black line) before tipping .The control parameter υ increases slowly with rate δ . In the ocean dynamics theforcing time scale is slow [23], so we set δ = 0.009 and ∆t = 0.05. The systemjumps to the lower branch from the upper branch, during which it goes througha smooth saddle node bifurcation. Because of the absolute value sign in the pro-cess (4.1), it can jump to another stable equilibrium on the lower branch, whichdoes not exist in the process (3.1) and (3.2).4.1.2 Estimate the linear decay rateAfter detrending, we are able to obtain the linear decay rate λ s that describes howfast y˜s relaxes back to the equilibrium after perturbed by noise. When the systemis stable, y˜s relaxes back quickly, while if υ is close to the bifurcation point, thesystem takes longer time to return, or even jumps away from the current state. The42above properties can be characterized through the linear decay rate λ sU and λ sACF .To obtain the linear decay rate λ sACF , we first calculate ACF coefficient αs byfitting it to a one-step autocorrelation function for y˜s in each sliding window. Thesize of a time window is again important and influences the result greatly. Sincethe 1-D model is different from the model in Chapter 3, we need to test differentwindow lengths again: we have tried four cases with w = N/2, w = N/3, w = N/4and w = N/5. According to the results of the linear decay rate, a sliding windowwith length w = N/4 works the best. For w = N/5 there are large fluctuations andwhen w = N/2,N/3, they are not able to capture the local information.After determining the window size w, suppose w = 2m+ 1, we measure αsjinside the jth window centered at υ j via the linear regressiony˜s,k+1 = αsj +σwky˜s,k (4.3)over k = j−m, ..., j +m using the Matlab package regress, which applies themethod of the least squares as described in Eq (3.9). After calculating αs in eachsliding window, we can work out λ sACF by λsACF =− logαs/∆t.As for λ sU from the effective potential well approximation, using the definitionof the potential well in Chapter 3, we have U = υy− y2 + 13 y3 for the systemgoing through the smooth saddle node bifurcation. The potential well Us for thedetrended data y˜s satisfies:dy˜s =−∂yUs(y˜s)dt+σdWtAnd the stationary Fokker-Planck equation (Eq (3.11)) ∂ p∂ t = L∗p= 0 derived fromdy˜s is12∂yy p(y˜s) = ∂y[−σ−2∂yUs(y˜s)p(y˜s)]. (4.4)After integration, it turns into12∂y p(y˜s) =−σ−2∂yUs(y˜s)p(y˜s)+σ−2c. (4.5)43We expand ∂yUs(y˜s) at the bottom of the well y˜s = 0, and in Chapter 3, we haveproved that after detrending ∂yUs(0) = 0, ∂yyUs(0) = λ sU , therefore:12∂y p(y˜s) =−λ sUσ−2y˜ · p(y˜s)+σ−2c (4.6)Next we can obtain the empirical density pemp(y˜s) in each sliding window withlength 2m+1 along the time series via the Matlab package ksdensity. For example,p j(y˜s) for the jth window centered at υ j is:p j(y˜s) =1(2m+1)hj+m∑i= j−mK(y˜s− y˜s,ih) (4.7)where h is default bandwidth h = σˆ( 43N )15 and σˆ is the standard deviation of thesample data, and K(·) is the kernel smoothing function [26]. We use the standardnormal distribution:K(x) =1√2pie−12 x2Then we compute the linear decay rate λ sU via nonlinear curve-fitting (Matlab pack-age lsqcurvefit). Because λ sUσ−2 is the product of two unknown variables, wedefine a new variable λ s′U and rewrite Eq (4.6) into12∂y p(y˜s) =−λ s′U y˜ · p(y˜s)+σ−2c (4.8)We compute λ s′U by solving the problem below using trust-region-reflective leastsquares algorithm (see Eq (3.16) in Chapter 3)minλ s′U‖ F(λ s′U , y˜s)−12∂y p(y˜s) ‖22= minλ s′U∑i(F(λ s′U , y˜s,i)−12∂y p(y˜s,i))2where F(λ s′U , y˜s,i) =−λ s′U y˜ · p(y˜s)+σ−2c. After obtaining λ s′U which equals λ sUσ−2,the empirical amplitude of noise can be estimated byσ2 =λ SACFλ s′U(4.9)Fig. 4.5 shows the linear decay rate λ sACF (blue) and λsU (red) of the time series44in Fig. 4.3. When tipping is pending, the linear decay rate decreases to zero. Thered curve and the blue curve have similar properties and both increase toward 1,indicating the system is only weakly stable.Figure 4.5: The linear decay rate of the detrended time series y˜s (in Fig. 4.4)with window length w = N/4: the red line is λ sACF and the blue line isλ sU .4.1.3 Nonlinearity information from the potential wellThe potential well of the process (4.1) when y > 0 isUs = υys− y2s +13y3sIt is a cubic function like Eq (3.18). So we use the same scheme that has been intro-duced in Section 3.1.3—approximate the potential well to a parabola. The potentialwell is related to the probability density via the Fokker-Planck equation (3.11), andthe stationary solution of it is:12∂yy p(ys) = ∂y[−σ−2∂yUs(ys)p(ys)]Gardiner [10] gives the solution of this equationUs =−σ22log p (4.10)Thus the empirical potential well can be estimated from the time series ys once wehave σ2 and p(ys). Recall λ s′U = σ−2λ sU , the noise amplitude can be estimated from45the linear decay rate λ sACF and λs′U . To estimate the empirical density p, we need tocompute it in each sliding window along the time series. That can be done with theMatlab package ksdensity: in the jth window, it solves (see Eq (3.14) for details)p j(ys) =1(2m+1)hj+m∑i= j−mK(ys− ys,ih) (4.11)using points ys, j−m,ys, j−m+1, ...,ys, j+m. The output pemp is a function of ys centeredat υ j and is evaluated at equally-spaced points ai that cover the range of the data inys. We name these points as (a1, p1),(a2, p2), ...,(aM, pM), then the empirical po-tential well centered at υ j can be determined by Eq (3.21), which is also a functionof ys, and we name these points as: (a1,U1), (a2,U2),...,(aM,UM).In Fig. 4.6, the gray surface is the potential well Us calculated from Eq (4.10)along the time window. The red line is a parabola fitting to the potential well fromthe last window. This parabola is computed by forming the Vandermonde matrixV and the parabola coefficients are equals to V/Us. The deviation of the approx-imated parabola from the empirical potential well is again measured by NRMSD(Eq (3.25)): for the time series in Fig. 4.3, it is 0.136, slightly larger than the stablelinear process 0.049 in Chapter 3.NRMSD =√∑Mi=1(Ui− (P1a21+P2a22+P3))2√M(Umax−Umin)The blue shadow below the potential surface is the time series ys and the minimumof the potential well gives the location of its equilibrium. If the process has alocally stable equilibrium or is a linear process, the red curve should overlap thepotential well boundary approximation based on Raymond [29]. The potential wellfor the current nonlinear process is rising above the parabola for the larger valueof ys while falling beneath the parabola for the smaller value of ys. According toEq (4.10), the part of the empirical potential well that falls below parabola indicatesthat there are more points in that direction, and a tipping to this direction is pending.46Figure 4.6: The nonlinear potential energy surfaces of the time series inFig. 4.3 using sliding window w=N/4. Us =−σ2emplog(pemp/2) whereσemp is estimated from Equation (3.22). The grey surface is the poten-tial energy, the red line is a parabola approximated to the potential wellvia Eq (3.24), the blue shadow below is the time series.4.1.4 SkewnessSkewness is again another indicator for tipping, and skewness γs is defined as:γs =1n ∑ni=1(ys,i− y¯s)(√1n ∑ni=1(ys,i− y¯s)2)3. (4.12)There is a Matlab package named skewness that implements this function tocompute γs. Fig (4.7) shows the skewness calculated from 100 realizations of pro-cess (4.1). Each realization is run from υ =−1 till υ = 0.9 or discarded if tippinghappens before υ = 0.9. Then for each discrete value of υ we have n = 100 differ-ent ys values coming from 100 realizations.From the figure, as υ moves closer to the bifurcation point, γs drops rapidlywhich indicates a larger likelihood of tipping. Besides, it is noticeable that γs(ys)always fluctuates for negative values and this indicates ys tends to jump in thenegative direction.47Figure 4.7: The estimated skewness of process (4.1) using 100 realizationsthat run from υ = −1 to υ = 0.9. Parameters: σ = 0.35, δ = 0.009,∆t = 0.05.4.2 Transition for the non-smooth bifurcationIn the previous section we have studied the system going through a smooth saddle-node bifurcation. But if υ decreases from a large value and undergoes a non-smooth bifurcation eventually, the dynamics of y are different.4.2.1 Time series and Gaussian kernel detrendingFirst we generate a time series (the black curve in Fig. 4.8) that starts at υ0 = 2with drifting rate δ = −0.009. The red line is the theoretical equilibrium, and thegreen line is the result of filtering the time series with a Gaussian kernel. Tippinghappens before η2 reaches the theoretical bifurcation point. So we pick out thetime series inside υ ∈ [0.2,2] and use yns to denote y in this time period. If onesupposes that the number of points in yns is N, then the Gaussian average Yns isYns(k∆t)k=1,2,...,N =∑Ni=1 Gk(i∆t)yns,i∑Ni=1 Gk(i∆t)whereGk(k=1,2,...,N)(t) =1√2pidexp(−12(t− k∆t)2d2)We use the same bandwidth d = 8 and ∆t = 0.05 as in the saddle node bifurcationcase. The right plot in Fig. 4.8 is the detrended time series y˜ns = yns−Yns.48Figure 4.8: [Left] The time series coming from the 1-D model Eq (4.1) withadditive noise. [Right] The detrended data y˜ns over Gaussian Kernel.Parameters: υ(0) = 2, σ = 0.35, δ = −0.009, ∆t = 0.05, bandwidthd = 8.4.2.2 Linear decay rateAfter detrending, we estimate the linear decay rate λ nsACF and λnsU . λ nsACF is computedfrom the autocorrelation function: in each window along the detrended data y˜ns, wefirst compute αns by linear regression:y˜ns,k+1 = αnsy˜ns,k +σwk.Calculation can be done with the Matlab package regress, in which it applies themethod of the least squares as shown in Eq (3.9). Then the linear decay rate λ nsACFof this process is related to αns via αns = exp(−λ nsACF∆t).For the linear decay rate λ nsU , we use the relationship between the detrendeddata y˜ns and its potential energy Uns. The potential well of the original process ynsis Uns = υyns− yns|yns|+ 13 y2ns|yns|. For the detrended time series, because we donot have an explicit form of ˙˜yns, we use its stochastic expressiondy˜ns =−∂yUns(y˜ns)+σdWtand derive the Fokker-Planck equation from this equation:12∂yy p(y˜ns) = ∂y[−σ−2∂yUns(y˜ns)p(y˜ns)],49which after integration, turns into12∂y p(y˜ns) =−σ−2∂yUns(y˜ns)p(y˜ns)+σ−2cSimilar to Section 3.1.2, λ nsU equals to ∂yyUns(0), so we expand ∂yU(y˜ns) at itsequilibrium—zero and get:12∂y p(y˜ns) =−σ−2y˜ns ·λ nsU · p(y˜ns)+σ−2cWe then combine σ−2λ nsU into one term, name it λ ns′U and the above function be-comes:12∂y p(y˜ns) =−y˜ns ·λ ns′U · p(y˜ns)+σ−2c (4.13)To obtain the empirical density p(y˜ns), we implement the kernel smoothing func-tion (Matlab’s ksdensity) in each window. For example, for the jth window withdata yns, j−m,yns, j−m+1, ...,yns, j+m centered at υ j,p j(y˜ns) =1(2m+1)hj+m∑i= j−mK(y˜ns− y˜ns,ih)where h is the bandwidth h = σˆ( 43N )15 , σˆ is the standard deviation of the sampledata and K(·) is the kernel smoothing function K(x) = 1√2pi e− 12 x2 . Once we havep(y˜ns), we use the nonlinear curve-fitting in the least-squares sense to get the coef-ficients λ ns′U and σ−2c by solving the problem as introduced in Eq (3.16)minc‖ F(λ ns′U , y˜)−12∂y pemp ‖22= mincj+m∑i= j−m(F(λ ns′U , y˜i)−12∂y pemp,i)2where F(λ ′U , y˜) =−σ−2y˜ns ·λ nsU · p(y˜ns)+σ−2cUnlike the saddle-node process in Fig. 4.5, the linear decay rate shown inFig. 4.9 does not have an obvious trend when approaching the non-smooth bi-furcation. λ nsACF and λnsU do not increase very much as υ decreases towards zero(non-smooth bifurcation point), and that means, the rate of the attraction for thestatic state does not vary much. Therefore, the linear decay rate is not able to50provide valuable information in this case.Figure 4.9: The linear decay rate of the detrended time series in Fig. 4.8:λ nsACF (red) from the autocorrelation function and λnsU (blue) from thepotential well. The window length for the linear and nonlinear analysisis w = N/4.4.2.3 Nonlinearity information from the potential wellThe theoretical potential well U of the time series (4.1) isUns = υyns− yns|yns|+ 13y2ns|yns|To obtain the empirical potential well of the time series in Fig. 4.3, we derive theFokker-Planck equation (3.11) of it:12∂yy p(yns) = ∂y[−σ−2∂yUns(yns)p(yns)],for which Gardiner [10] gave the stationary solution asUns =−σ22log p (4.14)Therefore, the potential well can be computed if σ2 and p are known. The noiseamplitude can be estimated through the linear decay rate λ nsACF and λns′U if one re-members that λ ns′U = σ−2λ nsU . To estimate the empirical density pemp, we obtain the51empirical density pemp in each sliding window via the probability density estima-tion. That can be done with the Matlab package ksdensity: in the jth window withpoints yns, j−m,yns, j−m+1, ...,yns, j+m, the package solves (see Eq (3.14) for details)pemp(yns) =1(2m+1)hj+m∑i= j−mK(yns− yns,ih) (4.15)The output pemp is a function of yns centered at υ j and is evaluated at equally-spaced points ai that cover the range of the data. We specify these points as(a1, p1),(a2, p2), ...,(aM, pM). Then from Eq (3.21), the potential well is also afunction of yns: (a1,U1), (a2,U2),...,(aM,UM). The nonlinearity is inferred fromthe difference between the potential well and the approximation parabola, which isestimated from the Vandermonde matrix: if one assumes the form of the potentialwell as P1a2+P2a+P3, P can be worked out as in Eq (3.24) bya21 a11 1a22 a12 1... ... ...a2M a1M 1P1P2P3=U1U2...UM (4.16)The deviation of the parabola from potential well is measured by NRMSD =√∑Mi=1(Ui−(P1a21+P2a22+P3))2√M(Umax−Umin) = 0.045.52Figure 4.10: The potential well of the time series in Fig. 4.8. The grey sur-face is the potential well, the red line is the approximated parabolicfit to the potential well in the last time window, and the blue shadowbelow them represents the time series yns. The potential energy Uns =−σ2emplog(pemp/2) where σemp is estimated from linear decay rate λ nsUand λ nsACF .Fig. 4.10 displays the potential well as a function of yns and υ from the time se-ries shown in Fig. 4.8. If we compare the potential well of the two time series withthe smooth and non-smooth bifurcation respectively, we have NRMSD = 0.045in the time series having non-smooth bifurcation (Eq (4.1)) and NRMSD = 0.136in the smooth counterpart (Eq (4.1)). NRMSD indicates the differences betweenthe approximated parabola and the empirical potential well, therefore indicates theprobability of tipping. One may conclude that the time series having non-smoothbifurcation is less likely to jump compared with the time series having smooth bi-furcation. But this is not true after a more careful consideration, which we willexplain subsequently.The time series having smooth bifurcation escapes at υ = 0.9 in Fig. 4.3. So werun 1000 realizations under the same conditions: σ = 0.35,δ = 0.009,∆t = 0.05.There are 380 out of 1000 realizations that have escaped by the time υ = 0.9. Wehence estimate the probability of tipping for the time series in Fig. 4.3 is 38%. Asfor the time series having non-smooth bifurcation in Fig. 4.8, it escapes at υ = 0.2.We again run 1000 realizations under same conditions, and there are 980 real-izations tipping before υ reaches 0.2, so the probability of tipping is 98%. The53non-smooth bifurcation actually has a higher probability of tipping. Therefore, theconclusion derived from NRMSD is incorrect.One possible reason that the NRMSD does not well predict the tipping is theabrupt change of behavior when crossing yns = 0. When this happens, the timeseries is actually a mixture of two different process: y˙ns = −υ + 2yns − y2ns +σWt (yns > 0) and y˙ns = −υ − 2yns + y2ns + σWt (yns < 0). This discontinu-ity influences the potential well as well. For the lower branch, we can not chooseone sign or the other because it crosses yns = 0. The potential well should beUns = υyns− yns|yns|+ 13y2ns|yns| (4.17)While for the time series that have the smooth bifurcation, we only need to considerthe upper branch, and so the potential well can be simplified asUs = υys− y2s +13y3s (4.18)The two potential wells differ from each other and lead to significantly differentbehaviors.Figure 4.11: The five potential curves of time series in Fig. 4.8 (left) and inFig. 4.3 (right). Red line—the theoretical potential, Eq (4.17) for theleft one and Eq (4.18) for the right one. Black line—the empiricalpotential curve, grey curve—the original time series generated fromEq (4.1), blue line—the theoretical equilibrium of the time series.54The empirical and the theoretical potential curves for the non-smooth (left) andthe smooth bifurcations (right) are shown in Fig. 4.11. The red line is the theoret-ical potential curve derived from Eq (4.17) and Eq (4.18). Below them, the greyshadow is the time series which follows theoretical equilibrium (blue line). Weonly draw the potential energy from 5 windows. For both empirical and theoreticalpotential curves, the minimum of the potential well gives the location of the equi-librium, where the red and black lines are coincided with each other. Furthermore,because both processes are quite stable at the very beginning of the time series, thedeviation is also small: NRSMD≈ 0.02 in both cases.The two potential wells in Fig. 4.11 are different. The right figure, in whichthe time series jumps to smaller ys eventually, has a more shallow potential wellon the larger side of ys (NRMSD = 0.136 in the final window), while for the leftplot, the potential well is deep at least near the equilibrium. The parabola approx-imation works until in the very end (NRMSD = 0.045 in the final window). Thisexplains why approximating the potential well by a parabola can not deduce theright conclusion in the 1-D model. The height of the potential well also explainsthe difference between the two types of bifurcations: in the right figure, the po-tential well for ys < yeq is very shallow, so the system easily escapes to anotherstate, while in the left figure, near yeq the well is still rather deep and then becomesshallow for yns larger than yeq, showing this system is likely to escape at the veryend.If we do not take the absolute value into effect, and generate the potential well:y˙ =−υ−2y+ y2+σWt (4.19)U = υy+ y2− 13 y3 (4.20)Fig. 4.12 shows the theoretical and empirical potential well for the time seriesy˙ns =−υ+2|yns|− yns|yns|+σWt (left), and y˙ =−υ−2y+ y2+σWt (right). Forthe right figure, its potential well is a cubic function and its potential well is thedeepest around equilibrium among these three cases—upper branch, lower branch55and lower branch without absolute value, indicating that the time series in the rightfigure is least likely to tip. In the right Fig. 4.12, the shape of the potential wellaround the equilibrium is parabolic and we can conclude that at the end of time,process 4.20 is still far away from its bifurcation. Hence the absolute value changesthe shape of potential well and causes discontinuity at zero.Figure 4.12: The theoretical and empirical potential well: the left is for thetime series y˙ns = −υ + 2|yns| − yns|yns|+ σWt and the right is fory˙ = −υ − 2y+ y2 +σWt . Red line—the theoretical potential curve,black line—the empirical potential curve, grey curve—the originaltime series, blue line—the theoretical equilibrium.4.2.4 SkewnessThe other method to detect nonlinearity is skewness. We generate 100 realizationsusing Eq (4.1) and for each realization, υ decreases from 2 to 0.5. If tipping hap-pens before υ = 0.5, we discard that realization. Hence at each discrete value of υ ,the skewness is calculated from 100 points coming from 100 different realizations.The calculation can be done with the Matlab package skewness:γns =1n ∑ni=1(yns,i− y¯ns)(√1n ∑ni=1(yns,i− y¯ns)2)3In Fig. 4.13, the skewness is positive all the time and increases greatly after υ = 1,indicating a tipping to the larger value of yns is pending.56Figure 4.13: The estimated skewness of yns using 100 realizations generatedby process (4.1) from υ = 2 to υ = 0.5. Parameters: υ(0) = 2, σ =0.35, δ =−0.009, ∆t = 0.05.Discontinuity at 0As υ moves closer to the non-smooth bifurcation, the absolute value comes intoeffect and changes the dynamics of system, and therefore changes the skewness. InFig. 4.14, the right column is obtained from the process y˙ = −υ−2y+ y2 +σWt ,while the left column is from y˙ns =−υ+2|yns|−yns|yns|+σWt . For all realizationsthat are used below, all yns fluctuate around the lower branch by the time υ = 0.3.The first row is the two time series generated from above two processes. Themiddle row shows the skewness of the two processes as a function of υ . Withoutabsolute sign, the skewness fluctuates around zero and indicates a stable state. Thelast row shows the probability density by picking out two realizations under twodifferent fixed υ . In the two plots, the left one with absolute sign always has alonger tail than the right one. And as υ decreasing, the difference becomes moreobvious which is consistent with the skewness: the difference between red and blueline when υ = 0.5 is larger than when υ = 0.4, showing the difference in skewnessbetween these two processes is getting larger. The whole Fig. 4.14 illustrates thediscontinuity at 0 and how it influences system skewness.57Figure 4.14: 1st row: the left figure is the time series y˙ns = −υ − 2|yns| −yns|yns|+σWt and the right figure is y˙ = −υ − 2y+ y2 +σWt . 2ndrow: the skewness of above process. 3rd row: the probability densityof the process generated at fixed υ : 0.4 and 0.5.4.3 ComparisonIn Chapter 4 we use the same methods that have been introduced in Chapter 3 toextract tipping information from the time series.The linear behavior is estimated from the linear decay rate λ . It indicates howfast the state relaxes back to the stable branch under a small disturbance. We intro-duce two different schemes based on the autocorrelation function and the potentialwell. λACF is obtained by fitting the detrended data to a one-step autocorrelationfunction Eq (4.3), and λ sU is derived by fitting the empirical probability densityto the Fokker-Planck Eq (4.6). Noting that λACF and λU measure the same thing,and the reason we calculate both of them is to get the empirical noise amplitudevia Eq (4.9). The difference between the Fig. 4.5—smooth bifurcation, and theFig. 4.9—non-smooth bifurcation, is that λ in Fig. 4.5 increase to one as υ ap-proaches the bifurcation point, while λ in Fig. 4.9 increases only slightly. So thatthe system in Fig. 4.5 is becoming weakly unstable, while the strength of the attrac-tion back to the equilibrium in Fig. 4.9 does not vary much. Therefore, the linear58decay rate gives the early warning signal in the smooth bifurcation case, while forthe lower branch it does not indicate tipping near the non-smooth bifurcation.Fig. 2.9 shows the equilibrium of the process on the lower branch without abso-lute sign (4.19). That explains why there is no discernable trend in the linear decayrate for the lower branch—when the parameter reaches the non-smooth bifurcationpoint in Eq (4.1), the dynamics in Fig. 2.9 is still far away from the smooth bifur-cation corresponding to the dynamics for y < 0.In addition to the linear behavior, the supplemental help from the nonlinearanalysis is also necessary. In Section 4.1.3, the tendency to jump is firstly mea-sured by the deviation of the approximating parabola (Eq (4.16)) from the empir-ical potential well (Eq (4.10)). The process going through the smooth bifurcationhas a larger deviation than the non-smooth counterpart, which indicates a largerlikelihood of tipping for the former. However, this conclusion is in contrast to ournumerical results. Recall that this method is based on Sieber [30], where it com-pares the likelihood of tipping between a quadratic and a cubic form potential wellvia approximating them to a parabola. However, as we’ve shown in Fig. 4.11, nowthe potential well of the lower branch has a quartic form, and it is flat only at thevery end. The shape of the potential well of the upper branch is not well approx-imated by a parabola away from equilibrium, because the potential well is gettingflat more gradually towards the negative value of ys. Fig. 4.12 further explains thediscontinuity at yns = 0 for the time series having the non-smooth bifurcation bycomparing the potential well between Eq (4.17) and Eq (4.20).Another way to measure nonlinearity is via the skewness (Eq (4.12)) whichdisplays a skewed data distribution. In Fig. 4.7 and Fig. 4.13, as υ approaches thebifurcation point, the magnitude of skewness increases which indicates a strongernonlinearity. Moreover, the sign of γ guides the direction that y tends to jump to:for the smooth bifurcation, γ < 0, thus y is likely to tip to a smaller value, while forthe non-smooth bifurcation, γ > 0, thus y is likely to tip to a larger value. For bothsmooth and non-smooth bifurcation, the variance of skewness near the bifurcationreveals a strong nonlinearity nearby. Furthermore, compared to the more gradual59change of γs, γns does not vary much until υ is close the the bifurcation point.60Chapter 5The 2D Stommel model withsmooth and non-smoothbifurcationsIn Chapter 3 and Chapter 4, we have tested the methods from Sieber [30]. In thischapter we are going to apply the same methods, including linear decay rate, po-tential well and skewness in the 2D Stommel model with additive noise. Comparedto the 1D system, the Stommel model is two dimensional and thus there are somedifficulties in applying the same methods. For example, we can easily derive thetheoretical form of the potential well in the 1D hysteresis model, and then comparethe shape of the potential well between the empirical one and the theoretical one.However, it is hard to calculate the theoretical potential well with an extra variablein the Stommel model. As for skewness, since T and ψ interact with each other,we have to consider fluctuations not only in T but also ψ . Hence there are somesignificant differences when applying the same approaches in the Stommel model.To better illustrate, we generate two time series with different noise amplitudes —one is set to be comparable with the 1D model (see Appendix(A.1) for calculation), and the other is much smaller than that. The system response ψ is determinedby two simultaneous processes of T and S, where Wt1 and Wt2 are two independentBrownian motions:61dT = (η1−T (1+ |T −S|))dt+σ1dWt1dS = (η2−S(η3+ |T −S|))dt+σ2dWt2dη2 = εdtψ = T −S.(5.1)5.1 Location of the smooth saddle node bifurcationIn this section, we apply the linear decay rate, potential well, and skewness toderive the early warning signals as the system approaches the smooth bifurcation.5.1.1 Time series and Gaussian kernel detrendingFirst we generate two time series by the process (5.1) using different sizes of noiseand then determine their trend via Gaussian kernel.Figure 5.1: The two time series coming from the process (5.1) with slowly in-creasing η2. For the left time series, the noise amplitudes σ1 =σ2 = 0.2.For the right time series, σ1 = σ2 = 0.04. The black curve is the timeseries, the red line is the theoretical equilibrium and the green line is theresult of filtering the time series with a Gaussian kernel. Parameters:η2(0) = 0, η1 = 3, η3 = 0.3, ε = 0.005, bandwidth d = 8.The black line in Fig. 5.1 is the time series generated within the process (5.1).We use ψs to represent the system varying through the smooth static bifurcation,and more specifically, ψs1 is the response with the larger noise amplitude σ1 =σ2 =0.2, which is comparable with the noise amplitude in the 1-D process, and ψs2 is62the one with the smaller noise σ1 = σ2 = 0.04. ψs1 sometimes crosses the unstablebranch near the bifurcation due to its large noise amplitude, while ψs2 fluctuatesaround the theoretical equilibrium. The drifting rate of the forcing parameter η2is ε = 0.005, which is also set to be comparable with 1-D model (A.1). The redline represents the static theoretical equilibrium. We select the time series beforetipping, and filter it with the Gaussian kernel to get the moving average:Ψs(k∆t)k=1,2,...,N =∑Ni=1 Gk(i∆t)ψs,i∑Ni=1 Gk(i∆t), k = 1,2, ...Nwhere N is the total number of points before tipping. Gk(·) is Gaussian kernel(Eq (3.3)) and ψs,i is the ith point of ψs. The detrended time series ψ˜s can then beestimated by ψ˜s = ψs−Ψs.5.1.2 Linear Decay RateAfter detrending, we are able to obtain the linear decay rate ΛsACF and ΛsU fromthe detrended data via fitting to the autocorrelation function and the potential wellrespectively.ΛsACF is computed by fitting it to a one-step autocorrelation function in eachsliding window along the detrended data ψ˜s. We test different window lengthsagain as in Section 3.1.2 and Section 4.1.2, and according to the the performancesof both temporal resolution and coefficient quality, we use window length w=N/3for the time series with noise amplitude σ = 0.2. Recall that we have shown in theSection 4.1.3 and Section 3.1.3, the potential well use the same length of windowas in the linear decay rate. For the time series with noise amplitude σ = 0.04, no-tice that inside a same window, the time series ψs2 with a smaller noise has a moreobvious trend than ψs1 with a larger noise. This trend influences the probabilitydensity, and then affects the potential well according to Eq (4.10). Therefore, weneed a smaller window for the time series with smaller noise to ensure that the po-tential well is not affected by this trend, and in this case, we use a window lengthN/10.63Suppose w = N/3 = 2m+ 1, in the jth window centered at η2, j, we have thelinear regressionψ˜s,k+1 = αsjψ˜s,k over k = j−m, ... j+m−1and thus αsj is calculated via minimizing (see Eq (3.9) for details)j+m−1∑k= j−m(ψ˜s,k+1−αsjψ˜s,k)2or use the Matlab package regress. After calculating αs in each window, ΛsACF canbe computed by ΛsACF =− logαs/∆t.As for ΛsU , the empirical potential well Us(ψ˜s) for the detrended data satisfies:dψ˜s =−∂ψUs(ψ˜s)dt+σψdWt (5.2)The stationary Fokker-Plank equation (see Eq (3.11)) of Eq (5.2) is12∂ψψ p(ψ˜s) = ∂ψ [−σ−2ψ ∂ψUs(ψ˜s)p(ψ˜s)],which after integration, is12∂ψ p(ψ˜s) =−σ−2ψ ∂ψUs(ψ˜s)p(ψ˜s)+σ−2c,To relate the above equation to ΛsU , we expand ∂ψUs(ψ˜s) around the equilibriumψ˜s = 0 of the detrended data. Since ∂ψUs(0˜s) = 0 and ∂ψψUs(0˜) = ΛsU ,12∂ψ p(ψ˜s) =−σ−2ψ ΛsU ψ˜ · p(ψ˜s)+σ−2c (5.3)We first estimate the empirical density p(ψ˜s) in each sliding window. In the jthwindow with length w = 2m+ 1, the empirical probability centered at η2, j is (or64using the Matlab package ksdensity)p j(ψ˜s) =1(2m+1)hj+m∑i= j−mK(ψ˜s− ψ˜s,ih)where h is the default bandwidth h = σˆ( 43N )15 , σˆ is the standard deviation of thesample data, and K(·) is the kernel smoothing function defined in Eq (3.14). Weuse Λs′U to represent σ−2ψ ΛsU . Then Λs′U can be computed by non-linear curve fit-ting using trust-region-reflective least squares algorithm or via the Matlab packagelsqcurvefit:minΛs′U‖ F(Λs′U , ψ˜s)−12∂ψ p(ψ˜s) ‖22= minΛs′U∑i(F(Λs′U , ψ˜s,i)−12∂ψ p(ψ˜s,i))2where F(·) is the right side of Eq (5.3).Unlike the model (4.1) and model (3.1), the system response ψs is influencedby both T and S. So that even though we have the noise of T and S in the artificialtime series, we still do not have σψ , which is needed to computed ΛsU . However,ΛsUσ−2ψ is not worthless, it will be used in the next section.Fig. 5.2 shows the results of the linear decay rate. The linear behavior of timeseries with different noise amplitude in Fig. 5.1 have similar properties. They bothincrease toward 1 as η2 approaches the smooth bifurcation point, indicating Λ isdecreasing to 0 and the system is only weakly stable.65Figure 5.2: The linear decay rate of process (5.1) with larger noise—(a), andsmaller noise—(b). Window length w = N/3 in the left one and w =N/10 in the right one.5.1.3 Nonlinearity information from the potential wellThe explicit form of the potential well Us(ψs) of the process (5.1) is difficult to get.However, we have the stochastic expression of ψs which isdψs =−∂ψUs(ψs)dt+σψdWtThe potential well is related to the probability density via the Fokker-Planck equa-tion (3.11), which has the stationary solution12∂ψψ p(ψs) = ∂ψ [−σ−2ψ ∂ψUs(ψs)p(ψs)]The solution of the above equation given by Gardiner [10] isUs =−σ2ψ2log p(ψs) (5.4)Thus the empirical potential well of ψs can be estimated once we have σ2ψ andp(ψs). We already know σ2ψ = ΛsACF/Λs′U in the Section 5.1.2. To estimate theempirical density p(ψs), we need to compute it in each sliding window along thetime series (or use the Matlab package ksdensity): in the jth window centered at66η2, j, we solve (see Eq (3.14) for details)p j(ψs) =1(2m+1)hj+m∑i= j−mK(ψs−ψs,ih)using points ψs, j−m,ψs, j−m+1, ...,ψs, j+m. Then the empirical potential well Us(ψs)centered at η2, j can be determined via Eq (5.4).The results of the potential wells are shown in Fig. 5.3. The left figure isgenerated from the time series with larger noise, while the right one has the smallernoise. The gray surface is the empirical potential well Us(ψs) calculated fromEq (5.4). The maximum of Us(ψs) in the right figure is much smaller than the leftone because of its much smaller noise amplitude. The blue curve below Us(ψs) isthe time series as shown in Fig. 5.1. The red line is the approximated parabola tothe potential well in the last time window. This parabola is computed by formingthe Vandermonde matrix V as shown in Eq (3.25). The approximation differenceNRMSD of the left plot is 0.216 whereas for the right figure it is 0.18, which leadsto the conclusion that the left time series (larger noise) is more likely to escape.Furthermore, the empirical potential well falls beneath the red line at the negativevalue side and indicates tipping is pending to that direction.Figure 5.3: The potential wells of the process (5.1). The left figure corre-sponds to the left time series in Fig. 5.1, and the right figure correspondsto the right time series in Fig. 5.1. Here Uemp =−σ2ψ log(pemp/2) whereσψ is estimated based on Equation (3.22). Sliding window w = N/3 inthe left figure, and w = N/10 in the right figure.675.1.4 SkewnessThe other method to study nonlinearity is skewness. For the two amplitudes ofnoise, we generate 500 realizations for each of them using Eq (5.1) and for eachrealization, η2 increases from 0 to 1.22. If tipping happens before η2 = 1.22,we discard that realization. Hence for each specific value of η2, the skewness iscalculated from 500 points that come from 500 different realizations. Suppose wehave n points (n=500 in this case), the skewness γs is (Matlab package skewness)γs =1n ∑ni=1(ψs,i− ψ¯s)(√1n ∑ni=1(ψs,i− ψ¯s)2)3Fig (5.4) shows the skewness of the time series (5.1). The left figure is with thenoise amplitude σ = 0.2, and the right figure is with σ = 0.04. The main differencebetween them is that γs of ψs1 decreases gradually and drops rapidly when η2 isclose to the bifurcation point, while γs of ψs2 does not vary much until very closeto the bifurcation point. This is because ψs2 does not have the tendency to jumpuntil very end as shown in Fig. 5.1. Overall, the variation of γs is a precursor oftipping and it works well in the branch with the smooth bifurcation, especially withthe larger noise.Figure 5.4: The skewness of the process (5.1) with slowly increasing η2. Theleft figure is generated with the larger noise and the right figure is gen-erated with the smaller noise.685.2 Location of the non-smooth bifurcationIn this section we test the same methods — linear decay rate, potential well andskewness in the non-smooth bifurcation case.5.2.1 Time series and Gaussian kernel detrendingFig. 5.5 shows the two time series coming from the process (5.1). The left figurehas noise amplitude σ1 = σ2 = 0.2, and the right figure has σ1 = σ2 = 0.04 as inSection 5.1, while for this time, η2 decreases from a larger value and passes throughthe static non-smooth bifurcation. In the left figure, sometimes ψns1 crosses zeronear the bifurcation due to its large noise amplitude, while in the right figure, ψns2does not cross 0 until η2 is very close to the bifurcation point.Figure 5.5: The two time series coming from the process (5.1) when η2 isslowly decreasing. The left figure has the noise amplitude σ1 = σ2 =0.2, and the right figure has σ1 = σ2 = 0.04. Parameters: η2(0) = 2,η1 = 3, η3 = 0.3, ε = 0.005, ∆t = 0.05.In Fig. 5.5, the black curve is the time series, the red line is the theoretical equilib-rium, and the green line is the Gaussian moving average Ψns:Ψns(k∆t)k=1,2,...,N =∑Ni=1 Gk(i∆t)ψns,i∑Ni=1 Gk(i∆t), whereGk(k=1,2,...,N)(t) =1√2pidexp(−12(t− k∆t)2d2)then the detrended time series is ψ˜ns = ψns−Ψns695.2.2 Estimate the linear decay rateThe linear decay rate ΛnsACF and ΛnsU of the the detrended time series ψ˜ns describeshow fast ψ˜ns relaxes back to the equilibrium 0 after perturbed by noise.ΛnsACF can be computed by fitting the ACF coefficient αns to the autocorrelationfunction in each sliding time window along y˜ns:ψ˜ns,k+1 = αnsψ˜ns,k +σψwk.The above calculation can be done with the Matlab package regress, in which itapplies the method of the least squares as shown in Eq (3.9). Then the linear decayrate ΛnsACF is related to αns via αns = exp(−ΛnsACF∆t).For the linear decay rate ΛnsU , we use the stochastic expression of the potentialwell of the detrended data ψ˜nsdψ˜ns =−∂ψUns(ψ˜ns)+σψdWtand derive the Fokker-Planck equation from this equation:12∂ψψ p(ψ˜ns) = ∂ψ [−σ−2ψ ∂ψUns(ψ˜ns)p(ψ˜ns)],which after integration, turns into12∂ψ p(ψ˜ns) =−σ−2ψ ∂ψUns(ψ˜ns)p(ψ˜ns)+σ−2cSimilar to the Section 3.1.2, ΛnsU equals to ∂ψψUns(0), so we expand ∂ψU(ψ˜ns)around the equilibrium ψ˜ns = 0 and get:12∂ψ p(ψ˜ns) =−σ−2ψ ψ˜ns ·ΛnsU · p(ψ˜ns)+σ−2ψ cWe then combine σ−2ψ ΛnsU into one term, name it Λns′U and the above function be-comes:12∂ψ p(ψ˜ns) =−ψ˜ns ·Λns′U · p(ψ˜ns)+σ−2ψ c (5.5)70To obtain the empirical density p(ψ˜ns), we implement the kernel smoothing func-tion (Matlab ksdensity) in each window. For example, in the jth window with dataψns, j−m,ψns, j−m+1, ...,ψns, j+m centered at η2, j,p j(ψ˜ns) =1(2m+1)hj+m∑i= j−mK(ψ˜ns− ψ˜ns,ih)where N is the number of points from the unknown distribution, h is the bandwidthh= σˆψ( 43N )15 , σˆψ is the standard deviation of the sample data and K(·) is the kernelsmoothing function K(x) = 1√2pi e− 12 x2 . Once we have p(ψ˜ns), we use the nonlinearcurve-fitting in the least-squares sense to get the coefficients Λns′U by solving theproblem as introduced in Eq (3.16)minc‖ F(Λns′U , ψ˜)−12∂ψ pemp ‖22= mincj+m∑i= j−m(F(Λns′U , ψ˜i)−12∂ψ pemp,i)2where F(Λ′U , ψ˜) = ψ˜ns ·Λns′U · p(ψ˜ns)+σ−2ψ cHowever, as we have pointed out in the Section 5.1.2, we are not able to com-pute ΛnsU because σψ is not known. Thus Fig. 5.6 only shows the results of thelinear decay rate ΛnsACF . (a) is the linear decay rate Λns1ACF of the process (5.1) withnoise amplitude σ = 0.2. (b) shows the linear decay rate Λns2ACF of the process (5.1)with noise amplitude σ = 0.04. Compared with the saddle node process that havea clear increasing trend over the range of η2 shown in Fig. 5.2, ΛnsACF of the lowerbranch does not have a gradual increase. The linear decay rate is not able to pro-vide much information at least in the first half period.We may also notice that the values of the linear decay rate in Fig. 5.6 is largerthan Fig. 5.2. As proved in Section 2.1.3, the system that fluctuates around theupper branch has two negative eigenvalues, then the solution should decay to theequilibrium exponentially. While for the system that fluctuates around the lowerbranch, it has two complex eigenvalues with negative real parts. Therefore, thesolution of the process on the lower branch is spiral back to the equilibrium, whichleads to a smaller decay rate Λns than Λs.71Figure 5.6: The linear decay rate of the detrended time series ψ˜ns. (a) is fromthe time series with the larger noise and (b) is from the smaller noise.5.2.3 Nonlinearity information from the potential wellThe stochastic expression of ψns related to the potential well Uns(ψns) isdψns =−∂ψUns(ψns)dt+σψdWtThen the Fokker-Planck equation of above is12∂ψψ p(ψns) = ∂ψ [−σ−2ψ ∂ψUns(ψns)p(ψns)]Here p(ψns) is the marginal probability density for ψns. The above Fokker-Planckequation has the stationary solution [10]:Uns =−σ2ψ2log p(ψns) (5.6)σ2ψ can be estimated by the results in the Section 5.2.2 — σ2ψ = ΛnsACF/Λns′U , andp(ψns) can be estimated by(Eq (3.14) for details): in the jth window along ψnswith points [ψns,i]k+mi=k−m inside, p(ψns) in that window isp j(ψns) =1(2m+1)hj+m∑i= j−mK(ψns−ψns,ih)Then the empirical potential well Uns(ψns) is computed by (5.6).72The results of the potential wells are shown in Fig. 5.7. It exhibits two potentialwells of the two time series as shown in Fig. 5.5 respectively. The time series in theleft one has the larger noise σ1 = σ2 = 0.2 and the right one has the smaller noiseσ1 = σ2 = 0.04. The deviation of the red approximated parabola (Eq (4.16)) fromthe potential well in the last window is measured by NRMSD. The results of bothcases are about 0.07. Compared with the NRMSD shown in the Fig. 5.3, which are0.216 and 0.18, the potential well of the time series on the lower branch fits theapproximated parabola much better and therefore, NRMSD is not able to indicatetipping near the non-smooth bifurcation.Figure 5.7: The potential well of the process (5.1). The grey surface is thepotential well, the red line is the approximated parabolic fit to the po-tential well in the last time window, and the blue shadow below themrepresents the time series ψns1 and ψns2 that show in Fig. 5.5.The explicit form of the potential well Uns(ψns) is hard to derive because itevolves a two dimensional process. However, this problem can be simplified if wefix T, eg. set T = Teq under each specific η2. Recall that:T˙ = η1−T (1+ |ψ|)→ Teq = η11+ |ψ|andψ˙ = (η1−η2−T +η3T )−ψ|ψ|−η3ψTherefore, when ψ > 073Teq =η11+ψU˙(ψ) = (−η1+η2)+ η1(1−η3)1+ψ +ψ2+η3ψ(5.7)Hence, the potential of ψ when ψ > 0 isU(ψ+) = (−η1+η2) ·ψ+η1(1−η3) ln(1+ψ)+ 13ψ3+12η3ψ2In a similar way, we can derive the potential of ψ when ψ < 0:U(ψ−) = (−η1+η2) ·ψ−η1(1−η3) ln(1−ψ)− 13ψ3+12η3ψ2 (5.8)Unlike the ψs on the upper branch that is always positive, we can not choose onesign or the other for ψns. In the left Fig. 5.8, we draw both the empirical (black)and theoretical (red) potential well of ψns using Eq (5.7) and (5.8). The lowestpoint, which indicates the equilibrium, of the red and black line are coincided witheach other. The last red line, which is at η2 = 1.1, is the average tipping time ofψns. There is no black line (potential well) drawn at that time because each blackline is computed by the data within a window. Suppose at η2 = η2,k, the data in-side the kth window is from ψns,k−m to ψns,k+m, hence the last black line is locatedat about half the window before η2 = 1.1. The right panel of Fig. 5.8 shows theempirical (black) and theoretical (red) potential well of ψs on the upper branch us-ing Eq (5.7). The final thick red line is at η2 = 1.26, which is the average tippingtime of ψs. Compared with the thick red line in the left figure, the right one ismore shallow around the equilibrium, which leads to a smaller NRMSD from theparabola approximation. However, the well around the equilibrium on the left sideis deeper, so that we would have a pretty good parabola fit even though the tippingis about to happen. And that explains why parabola approximation of Uns(ψ) doesnot work on the lower branch.74Figure 5.8: The theoretical (red) potential well and the empirical (black) po-tential well.If we do not take the absolute value into effect, the system would look likedT = (η1−T (1−ψ))dtdψ = (η1−η2+ψ2−T +(T −ψ)η3)dt(5.9)Fig. 5.8 compares the theoretical potential well of ψns (5.1) (red) and the theoret-ical potential well of ψ (5.9) (blue). The blue well has a perfect parabolic shapearound the equilibrium, indicating that the system is still far away from the bifur-cation and is stable. Hence we may conclude that the absolute value sign inside thesystem (5.1) causes the discontinuity at zero and leads to the non-smooth bifurca-tion.75Figure 5.9: The theoretical potential well of ψns of the deterministic versionof Eq (5.1) (red) and Eq (5.9)(blue).The potential well of TUnlike the 1D system, the Stommel model is a system that involves two pro-cesses (5.1). Since the NRMSD of ψns can not reveal too much information, inthis section we are going to discuss the potential well of T .Recall that the Stommel model isdTdt= η1−T (1+ |ψ|)dψdt= η1−η2−ψ|ψ|−T +(T −ψ)η3so that like in the potential well of ψns, we compute the theoretical potential wellof Tns via two steps — ψns > 0 and ψns < 0.— When ψns > 076ψ˙ns = 0→ (η1−η2−Tns+η3Tns)−ψ2ns−η3ψns = 0ψns =−η3±√η23 +4(η1−η2−Tns+η3Tns)2(take the positive sign)ψns =−η3+√η23 +4(η1−η2−Tns+η3Tns)2T˙ns = η1−Tns−TnsψUns(T ) =−η1Tns+(12 −η34)T 2+5bTns(√a+bTns)3−2(√a+bTns)515b2where a = η23 +4η1−4η2,and b = 4(η3−1)— When ψns < 0ψ˙ns = 0→ (η1−η2−Tns+η3Tns)+ψ2ns−η3ψns = 0ψns =−η3±√η23 −4(η1−η2−Tns+η3Tns)2(take the negative sign)ψns =−η3−√η23 −4(η1−η2−Tns+η3Tns)2T˙ns = η1−Tns+TnsψnsUns(T ) =−η1Tns+(12 −η34)T 2+5bTns(√a+bTns)3−2(√a+bTns)515b2where a = η23 +4η2−4η1,and b =−4(η3−1)Fig. 5.10 shows the theoretical (left) and empirical (right) potential well of Tns.The NRMSD between the potential well in the last window and the approximatedparabola is NRMSD = 0.18, which is much larger than the NRMSD derived fromUns(ψ). One may also notice from the left figure, that as η2 approaches the non-smooth bifurcation, because ψns fluctuates up and below ψ = 0, if Tns gets insidethe region where the blue line and the red line overlaps, the shape of the empiricalpotential well is hard to be approximated by a parabola.77Figure 5.10: The left panel is the theoretical potential well of Tns. The redline is when ψns < 0, and the blue line is when ψns > 0. The rightpanel is the empirical potential well of Tns (grey) and the approximat-ing parabola (red).5.2.4 SkewnessFig (5.11) shows the distribution of the skewness γns of the time series generatedby Eq (5.1) with a negative ε viaγns =1n ∑ni=1(ψns,i− ψ¯ns)(√1n ∑ni=1(ψns,i− ψ¯ns)2)3where n is the number of the realizations we use, and n = 500 in this case. The toppanel is generated with noise amplitude σ = 0.2, and the bottom panel is generatedwith noise amplitude σ = 0.04. As shown in Fig. 5.5, the right time series whichhas noise amplitude σ = 0.04 fluctuates around the equilibrium quite steadily be-fore the bifurcation point. Hence γns of ψns1 increases earlier than γns of ψns2,and also have a larger variation than ψns2. Compared with the skewness γs in theSection 5.1.4, γns does not vary much until η2 is very close to the bifurcation point.78Figure 5.11: The skewness of ψns of the process (5.1) with slowly decreasingη2. The top panel is generated with the larger noise and the bottompanel is generated with the smaller noise.The skewness of TUnlike the 1D model in the Chapter 3 and Chapter 4, the Stommel model is a 2Dmodel and therefore the skewness of T is also worth being computed. Fig. 5.12displays the skewness of ψns1 and Tns1 respectively. They are computed using thesame 500 realizations and it is noticeable that the range of the skewness of Tns1 islarger than ψns1.79Figure 5.12: Comparison between the skewness of ψns1 and Tns1.Fig. 5.13 is the trajectories in the T and ψ direction when η2 = 1, and thered circle is the starting point which then follows the red line and moves backto the equilibrium (2.88,−0.05). The trajectory in the T direction has a largervariance than in the ψ direction which explains why the skewness of Tns has amore discernible trend than ψns.Figure 5.13: Trajectories in the direction of ψ and T on the lower branch atη2 = 1.80Chapter 6Prediction of non-smoothbifurcation via coherenceresonanceApart from the methods from Sieber [30], recall that in Section 2.1.3, we show thatthe eigenvalues of the Stommel model on the lower branch are two complex num-bers with negative real parts. As η2 approaches the non-smooth static bifurcationpoint, these eigenvalues also vary gradually, which has an impact on the way thesystem oscillates as well. Fig. 6.1 displays the oscillations of ψns generated bydTdt= η1−T (1+ |T −S|)dSdt= η2−S(η3+ |T −S|)ψ = T −Swith different fixed η2 and additive noise. The panels on the right show the cor-responding power spectral density. The top realization with η2 = 1.3 has a morediscernable periodic cycle than the realization below, and the top one has a clearerpeak in its power spectral density.Suppose the eigenvalues of Eq (1.9) are two complex numbers λ = −a± bi81(a,b > 0), then the solutions of ψ and T look like:[ψT]=C1e−at[cosbtsinbt]+C2e−at[sinbt−cosbt]where C1 and C2 are two constants. Thus the larger the magnitude of a, the fasterthe solution decays. We compute the eigenvalues of both cases via Eq (2.4), andhave a = 0.9 in the figure above and a = 3 in the figure below. With a = 3, ψnsdecays back to the equilibrium relatively rapidly and after that its oscillations aremainly dominated by noise. Thus we can not observe a clear periodic cycle.Figure 6.1: The oscillation of ψns generated with fixed η2 and its correspond-ing power spectral density. The figure above is generated with η2 = 1.3,which has eigenvalues −0.9± 1.3i , and the figure below is generatedwith η2 = 5, which has eigenvalues −3.0± 0.8i. Other parameters:η1 = 3, η3 = 0.3, σ = 0.05, ∆t = 0.05, tend = 30.With respect to the periodic cycle, Pikovsky [28] introduced a scheme thatcan detect the character of the underlying stochastic dynamics based on the powerspectral density. The power spectral density interprets the distribution of power into82the frequency components composing that signal. If the power is concentrated neara single frequency, the signal is said to be coherent. Coherence resonance describesthe coherent behavior such as rhythmic oscillations. The coherence resonance β ofa dynamics is defined as [9]β = hp(∆w/wp)−1, (6.1)where hp and ∆w are the height and the width of the spectrum peak at frequencywp of the power spectral density.6.1 The coherence resonance in different realizationsAs η2 approaches the non-smooth static bifurcation point, the real part a of theeigenvalue of (2.3) decreases in magnitude gradually (Eq (2.5)), so that the damp-ing of the oscillations is weaker. This indicates that the periodic cycle shouldbe observed more clearly and the power spectral density should have an obviouspeak. However, when η2 is close to the non-smooth static bifurcation, becauseψns crosses zero from time to time, the eigenvalues of the system are sometimescomplex numbers and sometimes real numbers. The trajectory then is affected bythe discontinuous dynamics, which may lead to an unclear peak in the power spec-tral density (PSD). Later, we are going to test three different cases to see if thecoherence resonance can give an early warning signal in each of them.6.1.1 η1 = 3,η3 = 0.3In the first example, we set η1 = 3 and η3 = 0.3 as in Chapter 5. The theoreticalnon-smooth static bifurcation point is at η2 = 0.9. We generate 300 realizationsat each discrete value of η2 ∈ [1,1.6], and in order to detect the cycle period, eachrealization is set to be run until t = 50, and we discard the realizations that haveescaped by t = 50. Since the system is very likely to jump near the bifurcationpoint, we use noise amplitude σ = 0.05.Fig. 6.2(b) shows the measure of coherence resonance β as a function of η2which is computed by Eq (6.1). Firstly we calculate the power spectral density83for each of the realizations. Then we compute the coherence resonance of eachrealization by Eq (6.1) after determining hp and ∆w as described below. Supposewe have N points in a realization (t1,ψns,1),(t2,ψns,2), ...,(tN ,ψns,N), we first needto compute its discrete Fourier transform. In this thesis, we use the Matlab packagepwelch which follows the below steps:Let vk, k ∈ R be the Fourier transform of this realization,vk =N∑j=1(ψns, j− ψ¯ns)(exp(−2pii/N))(k−1)( j−1)as k goes to infinity, the Fourier transform decreases to zero. After that, the powerspectral density (PSD) of that realization as a function of frequency k can be com-puted via [13]:PSD =∆tN· |vk|2Figure 6.2: (a) The bifurcation diagram of ψ when η1 = 3 and η3 = 0.3. (b)The measure of coherence resonance β of the stochastic dynamics (1.9).Other parameters: tend = 50, ∆t = 0.01, noise amplitude σ = 0.05.Suppose the highest point of that power spectral density is (wp,hp), that is,after computing PSD, we find the maximum PSD value and set it to be hp, andset its corresponding frequency to be wp. The width ∆w of the peak is defined asfollowed. Firstly we find two points (a,10−6) and (b,10−6), which means whenfrequency equals a or b, the power spectral density is 10−6. Then the width is de-fined as the difference between a and b. After having hp, wp, and ∆w, β can becomputed by Eq (6.1). In Fig. 6.2, (b) is the averaged coherence resonance over the84300 realizations, and it reaches the maximum value at η2 = 1.067, which has a rel-ative error of 18% compared with the non-smooth static bifurcation point η2 = 0.9.Figure 6.3: The probability densities of ψns computed using 500 realizationsat (a) η2 = 1.2, (b) η2 = 1.067, (c) η2 = 1. Other parameters: η1 = 3,η3 = 0.3.To understand the behavior of ψns around η2 = 1.067, Fig. 6.3 compares theprobability densities (PDF) of ψns at different values of η2. (a) is the PDF of ψnswhen η2 is far from the bifurcation. All the realizations at η2 = 1.2 are stable andfluctuate around the equilibrium on the lower branch. (b) shows the PDF of ψnswhen β reaches its highest value. By that time we can notice that a few realizationshave already made the transition. When η2 decreases to 1, most of the realizationsfluctuate around the equilibrium on the upper branch as shown in (c). Therefore,the maximum value of β indicates the time when ψns starts to make the transition.85Figure 6.4: The power spectral density of ψns at (a) η2 = 1.2, (b) η2 = 1.067,(c) η2 = 1.Then to further explain the trend of β , Fig. 6.4 shows the power spectral densityat different value of η2. We compute the the power spectral density at η2 = 1.2,η2 = 1.067 and η2 = 1 respectively. In (a), because η2 is far from the non-smoothstatic bifurcation point, its power can be concentrated around a single frequency.Thus its PSD has a clear single peak. When η2 = 1, there are two discernablepeaks in the power spectral density — one peak is at Frequency≈ 0 and the otherpeak is at Frequency ≈ 2. Recall that we need to subtract the equilibrium valueof ψns from the lower branch so that only the oscillations are analyzed. However,there are some irregular fluctuations near the non-smooth static bifurcation point,which create a combination of different frequencies in the power spectral density.As for (b), because η2 = 1.067 is when the system starts to make the transition,the peak at Frequency ≈ 0 is not as obvious as in (c). Moreover, compared with(b), the height of PSD in (c) is smaller, and the width is larger. Thus the coherenceresonance (6.1) decreases after η2 = 1.067. These figures in (6.4) together explainthe trend of β in Fig. 6.2(b), and indicate that as η2 approaches the non-smoothstatic bifurcation, the system starts to experience larger irregular fluctuations aboutthe lower equilibrium states.866.1.2 η1 = 6,η3 = 0.4Figure 6.5: (a) The bifurcation diagram when η1 = 6 and η3 = 0.4. (b) Themeasure of coherence resonance β computed by Eq (6.1). Other param-eters: tend = 50, ∆t = 0.01, noise amplitude σ = 0.05.In the second case, we set η1 = 6, η3 = 0.4, and the non-smooth bifurcation is atη2 = 2.4. As in the previous case, we also generate 300 realizations at each discretevalue of η2 ∈ [2.6,3]. Fig 6.5(b) shows the averaged coherence resonance β over300 realizations as a function of η2 and β reaches the highest value at η2 = 2.72.Compared with the theoretical bifurcation at η2 = 2.4, our prediction has a relativeerror of 13%.Figure 6.6: The probability densities of ψns computed using 500 realizationsat (a) η2 = 3, (b) η2 = 2.72, (c) η2 = 2.55. Other parameters: η1 = 6,η3 = 0.4, tend = 50.87Fig. 6.6 shows the probability densities (PDF) of ψns at different values of η2.(a) is the PDF of ψns when η2 is far from the bifurcation point, and all ψns fluctuatearound the equilibrium on the lower branch. (b) shows the PDF when β reachesits highest value. Like in Fig. 6.3, we see that a few realizations have made thetransition to the upper branch. By the time η2 decreases to 2.55, most of the real-izations escape from the lower branch and fluctuate around the equilibrium on theupper branch as shown in (c). Therefore, the maximum value of β in this case alsoindicates the time when ψns starts to make the transition.Figure 6.7: The power spectral density of ψns at (a) η2 = 3, (b) η2 = 2.72, (c)η2 = 2.55.To help understand the trend of β , Fig. 6.7 computes the power spectral densi-ties of ψns at η2 = 3, η2 = 2.72, η2 = 2.55 respectively. And similar to the resultsin Fig. 6.4, because η2 = 3 is far from the non-smooth static bifurcation point,there is a one single peak in (a). In the next figure, there is a hill rising abovefrom the left side (Frequency≈ 0) because ψns starts to make the transition by thetime η2 = 2.72. When η2 = 2.55, the irregular fluctuations near the non-smoothstatic bifurcation create a combination of different frequencies. There are two dis-cernable peaks in the power spectral density — one peak is at Frequency ≈ 0 andthe other peak is at Frequency ≈ 1.5. These figures further explain the trend of βin Fig. 6.5(b). The decrease of β after η2 = 2.72 indicates that as η2 approaches88the non-smooth static bifurcation, the system starts to experience larger irregularfluctuations about the lower equilibrium states.6.1.3 η1 = 8,η3 = 0.3Figure 6.8: (a) The bifurcation diagram when η1 = 8 and η3 = 0.3. (b) Themeasure of coherence resonance β of ψns of the process (1.9) computedby Eq (6.1). Other parameters: tend = 50, ∆t = 0.01, noise amplitudeσ = 0.05.In the last example, we set η1 = 8, η3 = 0.3, and the non-smooth bifurcation is atη2 = 2.4. Fig 6.8(b) shows the coherence resonance β as a function of η2. Whenη2 = 2.84, β reaches its highest value and then drops off. The relative error of ourwarning prediction is 18%.89Figure 6.9: The probability densities (PDF) of ψns computed using 500 re-alizations at (a) η2 = 3, (b) η2 = 2.84, (c) η2 = 2.7. Other parameters:η1 = 8, η3 = 0.3, tend = 50.To understand the behavior of ψns around η2 = 2.84, we again calculate theprobability density distribution of ψns using 500 realizations at different η2. Allthe realizations in Fig. 6.9(a) fluctuate around the lower branch, while a few ofthem escape to the upper branch when η2 = 2.84. By the time η2 = 2.7, because itis close to the bifurcation point, most of the realizations have made the transitions.Therefore the decrease of β is an early warning signal.90Figure 6.10: The power spectral densities of ψns computed using 500 real-izations at (a) η2 = 3. (b) η2 = 2.84. (c) η2 = 2.7. Other parameters:η1 = 8, η3 = 0.3, tend = 50.To understand the behavior of β in Fig. 6.8, Fig. 6.10 shows the power spectraldensity (PSD) of ψns at three different values of η2. (a) because η2 is far from thenon-smooth static bifurcation point, its power can be concentrated near a singlefrequency. When η2 = 2.84, the PSD looks similar to (a), but because there area few realizations escaping to the upper branch, we may notice that the PSD atFrequency ≈ 0 becomes larger than (a). Then when η2 = 2.7, there are two dis-cernable peaks in the power spectral density — one peak is at Frequency ≈ 0 andthe other peak is at Frequency ≈ 2. Recall that we need to subtract the equilib-rium value of ψns from the lower branch so that only the oscillations are analyzed.However, there are some irregular fluctuations near the non-smooth static bifurca-tion, which create a combination of different frequencies. These figures in (6.10)explain the trend of β in Fig. 6.8(b). To sum up, in the above three examples,the decrease of β gives us the warning signal for capturing the route to tippinggenerated through coherence resonance.6.2 Coherence resonance for smaller noiseFig. 6.11 shows the coherence resonance of ψns with a smaller noise amplitudeσ = 0.01. We take the case η1 = 8,η2 = 0.3 as an example. Compared with91Fig. 6.8(b), since the dynamics has a stronger coherence resonance behavior in thelarger noise case, the magnitude of β in Fig. 6.8(b) is larger than in Fig. 6.11. Be-sides, β reaches its maximum value at η2 = 2.5 with a relative error of 4% in thesmaller noise amplitude.Then we compare the relative distance between the average tipping time andthe time when β reaches its maximum value. With the smaller noise amplitude,the average tipping time is when η2 = 2.48 and β reaches its maximum valuewhen η2 = 2.5. With the larger noise amplitude, the average tipping time is whenη2 = 2.74 and β reaches its maximum value when η2 = 2.84. The relative distanceof the larger noise case is longer. Therefore, a larger noise is likely to generate alarger advance of the tipping point relative to the non-smooth bifurcation.Figure 6.11: The coherence resonance β as a function of η2 when η1 =8,η3 = 0.3, the noise amplitude σ = 0.01.In Chapter 6, we detect the early warning signal by the change of the coherenceresonance. This method is only applicable for a two or higher dimensional systembecause only a square system matrix can have complex eigenvalues which lead toa periodic cycle in the system response [9]. In Chapter 5, the red trajectories inFig. 5.13 have spiral shapes due to the complex eigenvalues as calculated in theSection 2.1.3. As η2 approaches the non-smooth static bifurcation, the real part ofthe complex eigenvalues decreases in magnitude, and therefore change the powerspectral density. This change together with the shift of the system equation (1.9)92between the positive and negative ψns gives a varying coherence resonance, andaccording to the examples that we test, gives us the early warning signal.93Chapter 7Conclusion7.1 Summary of resultsThe main task of this thesis is to identify practical early warning signals that canbe used for both the smooth and non-smooth transitions corresponding to tippingnear the bifurcations in the Stommel model. The Stommel model depicts the over-turning in the Atlantic circulation and after being added with noise, it can be usedto simulate ocean climate in the Atlantic circulation. The Stommel model goesthrough different kinds of bifurcations when the input parameter η2 increases ordecreases. In the Appendix (A.1), we classify the smooth bifurcation in the Stom-mel model as a saddle node bifurcation via analytical calculation, and define asimplified 1D hysteresis model (2.8) that has similar properties as the Stommelmodel.Since the saddle node bifurcation has been well studied, we first review the lin-ear and nonlinear analysis—linear decay rate, potential well and skewness [30] inthe saddle node normal form in chapter 3. All of the methods above have been stud-ied in different settings to predict the saddle-node bifurcation. Firstly, a decrease ofthe linear decay rate to zero indicates a tipping instability [14] [19]. In terms of thepotential energy, in the immediate neighborhood of a stable equilibrium, the poten-tial energy can be fit quite well by a parabola. Hence when the system is stable andfar away from the bifurcation point, its potential energy should be approximated94by a parabola. But the parabolic shape of it will become increasingly perturbedby nonlinear features. Therefore we may use the deviation of the potential energyfrom its approximated parabola to indicate the likelihood of tipping. As for theskewness of the empirical density of the time series, its trend is likely to be influ-enced by the closeness of a bifurcation. Therefore, the results of these methodsare able to give us early warning signals when analyzing the smooth saddle nodebifurcation.We then apply these methods in the 1D hysteresis model to predict the tippingin chapter 4. The results of the smooth bifurcation case and the non-smooth bifur-cation case are widely different. The linear decay rate has an identifiable trend tozero in the smooth bifurcation case. While in the non-smooth bifurcation case, itdoes not have a discernable variation. As for the potential well, for the smooth bi-furcation case, the large deviation of the approximated parabola from the empiricalpotential well indicates that ψs has a tendency to escape as η2 moves closer to thesmooth bifurcation point. Whereas for the non-smooth case the deviation is muchsmaller and can not show any nonlinear features. Regarding the skewness, un-like the gradual change in the smooth bifurcation case, the skewness moving alongthe lower branch does not vary much until the parameter is very close to the non-smooth bifurcation point. Based on above analysis, the methods from Sieber [30]can give us the early warning signal when the system approaches the smooth bifur-cation, but does not work well for the non-smooth bifurcation.Next in chapter 5, we apply the same methods to the 2D Stommel model. Sim-ilar to the results in the 1D hysteresis model, those methods work well for thesmooth bifurcation, whereas for the non-smooth bifurcation, they can not giveenough information. However, unlike the 1D hysteresis model, the dynamics ofthe variable T is also worth evaluating. The trajectories in the Fig. 5.13 shows thatwhen the system moves back to the equilibrium on the lower branch, the path inthe T direction has a larger variation than in the ψ direction. Hence in the Sec-tion 5.1.3 and (5.1.4), we compute the potential well and the skewness of Tns, andfind out that the results of both methods have a more obvious warning signal thanψns.95Another important property of the dynamics on the lower branch of the 2Dmodel is the coherence resonance. Because the square system matrix when ψ < 0has two complex eigenvalues with negative real parts, there is the potential fornoise to sustain these oscillations with frequency near the imaginary part of theeigenvalues. Hence in chapter 6 we derive an early warning signal by consideringthe variation in the coherence resonance. Our study of the non-smooth bifurcationsuggests that as the system moves towards the non-smooth bifurcation, the coher-ence resonance will first increase and then decrease. The start of the decrease isan early warning signal that the system becomes unstable and has the tendency totip. First, the power spectral density is varying because of the change in the com-plex eigenvalues. Then as η2 moves closer to the bifurcation, because ψns crosseszero from time to time, the eigenvalues of the system are shifting between complexnumbers and real numbers. These two factors together give us the early warningsignal and make up for the insufficient warning information for the non-smoothtipping.7.2 Future directionsThere are some questions that may need to be answered in future work. For exam-ple, we might be able to study the variation of the behavior of the potential wellat different times over the time series. Right now the probability of tipping can bemeasured by the deviation of the potential well from its approximated parabola.And we only compute the deviation in the last time window. Then we comparedifferent deviations: in Section 3.1.3, we compare between the saddle node bi-furcation model and the linear model; in Chapter 4 and Chapter 5, we comparebetween the smooth bifurcation case and the non-smooth bifurcation case. How-ever, if we do not have another comparative time series, it also makes sense to dothe comparison between the earlier and later time within one process. One possi-ble concrete method is that we set an threshold, and once the change of NRMSDexceeds this value, the time series is very likely to escape from its current state.96Another possible improvement of this model is readjusting the parameters byreal data. In the Stommel model, we estimate the values of η1,η2,η3 by Table (2.1)roughly, and set the varying rate of the input parameter referring to [30]. In thepast 50 kyr, there are at least five abrupt flippings of the thermohaline circulationsystem [1]. Therefore, it might be worthwhile testing the usefulness of this modelby using the real MOC data. If the original setting is unsatisfactory, for example,tuning is one of the intrinsic and fundamental methods to optimize parameter [16].It chooses parameter values by minimizing the deviation of the model output fromhistorical data, or reducing it to an acceptable rage. By exploring a variety of sizesof these parameters, we can achieve the most realistic numerical model.In this thesis we focus on the Stommel model because this model can be suc-cessfully used to understand the effects of the freshwater flux and heat flux on theNorth Atlantic ocean circulation. However, there is no universal model that can beapplied to answer all the questions when tackling the ocean climate. Right nowthere are many other models that have been widely discussed, such as the 3-boxmodel by Rooth [20] and energy balance models (EBMs) by Budyko [34]. Thedifferences in their forms may lead to different simulation results, and thus it isalso worth experimenting these methods developed in this thesis in other models.For example, when we deal with a high dimensional model (2D or 3D), the firstmove is to determine the type of the transition of the variable that we study. Inthe Stommel model, we are able to conclude that the smooth bifurcation of ψ is asaddle node bifurcation, but for a system that involves more variables, sometimesit is hard to simplify the form and decide which variable shall be analyzed. Besidesit is also hard to derive the theoretical form of the potential well in a high dimen-sional system when studying the nonlinear properties. For other methods that donot need explicit forms, Lenton [17] investigated the thermohaline circulation viaanalyzing the linear decay rate of the three models — CLIMBER2, GENIE1, andGENIE2. Of these three models, the result of the most complex GENIE2 givesthe most clear signal of tipping. Accordingly, there are many aspects in a model,such as the time scales and the ocean components taken into account, affect theperformance of an approach, and the approaches developed in this thesis should beapplied in other models after considering these differences.97Bibliography[1] Modeling the thermohaline circulation.http://www3.geosc.psu.edu/∼dmb53/DaveSTELLA/THC%20model.pdf. →page 97[2] P. Ashwin, S. Wieczorek, R. Vitolo, and P. Cox. Tipping points in opensystems: bifurcation, noise-induced and rate-dependent examples in theclimate system. Phil. Trans. R. Soc. A, 370(1962):1166–1184, 2012. →page 3[3] P. Ashwin, S. Wieczorek, R. Vitolo, and P. Cox. Tipping points in opensystems: bifurcation, noise-induced and rate-dependent examples in theclimate system. Phil. Trans. R. Soc. A, 370(1962):1166–1184, 2012. →page 17[4] A. W. Bowman and A. Azzalini. Applied smoothing techniques for dataanalysis: the kernel approach with S-Plus illustrations, volume 18. OUPOxford, 1997. → page 29[5] G. E. Box, G. M. Jenkins, G. C. Reinsel, and G. M. Ljung. Time seriesanalysis: forecasting and control. John Wiley & Sons, 2015. → page 25[6] T. F. Coleman and Y. Li. An interior trust region approach for nonlinearminimization subject to bounds. SIAM Journal on optimization, 6(2):418–445, 1996. → page 29[7] T. F. Coleman and Y. Li. An interior trust region approach for nonlinearminimization subject to bounds. SIAM Journal on optimization, 6(2):418–445, 1996. → page 29[8] H. A. Dijkstra. Nonlinear climate dynamics. Cambridge University Press,2013. → pages 6, 1398[9] H. Gang, T. Ditzinger, C. Ning, and H. Haken. Stochastic resonance withoutexternal periodic force. Physical Review Letters, 71(6):807, 1993. → pages83, 92[10] C. W. Gardiner. Handbook of stochastic methods for physics, chemistry andthe natural sciences, p 415, 2004. → pages 33, 45, 51, 66, 72[11] V. Guttal and C. Jayaprakash. Changing skewness: an early warning signalof regime shifts in ecosystems. Ecology letters, 11(5):450–460, 2008. →page 37[12] V. Guttal and C. Jayaprakash. Spatial variance and spatial skewness: leadingindicators of regime shifts in spatial ecological systems. TheoreticalEcology, 2(1):3–12, 2009. → page 37[13] G. Heinzel, A. Ru¨diger, and R. Schilling. Spectrum and spectral densityestimation by the discrete fourier transform (dft), including a comprehensivelist of window functions and some new at-top windows. 2002. → page 84[14] H. Held and T. Kleinen. Detection of climate system bifurcations bydegenerate fingerprinting. Geophysical Research Letters, 31(23), 2004. →pages 24, 94[15] H. Held and T. Kleinen. Detection of climate system bifurcations bydegenerate fingerprinting. Geophysical Research Letters, 31(23), 2004. →page 25[16] F. Hourdin, T. Mauritsen, A. Gettelman, J.-C. Golaz, V. Balaji, Q. Duan,D. Folini, D. Ji, D. Klocke, Y. Qian, et al. The art and science of climatemodel tuning. Bulletin of the American Meteorological Society, 98(3):589–602, 2017. → page 97[17] T. Lenton, V. Livina, V. Dakos, E. Van Nes, and M. Scheffer. Early warningof climate tipping points from critical slowing down: comparing methods toimprove robustness. Phil. Trans. R. Soc. A, 370(1962):1185–1204, 2012. →pages 23, 97[18] T. M. Lenton, H. Held, E. Kriegler, J. W. Hall, W. Lucht, S. Rahmstorf, andH. J. Schellnhuber. Tipping elements in the earth’s climate system.Proceedings of the national Academy of Sciences, 105(6):1786–1793, 2008.→ page 399[19] V. N. Livina and T. M. Lenton. A modified method for detecting incipientbifurcations in a dynamical system. Geophysical Research Letters, 34(3),2007. → page 94[20] V. Lucarini and P. H. Stone. Thermohaline circulation stability: A boxmodel study. part i: Uncoupled model. Journal of climate, 18(4):501–513,2005. → page 97[21] V. Lucarini and P. H. Stone. Thermohaline circulation stability: A boxmodel study. part i: Uncoupled model. Journal of climate, 18(4):501–513,2005. → page 16[22] I. D. Mayergoyz. Mathematical models of hysteresis and their applications.Academic Press, 2003. → page 4[23] J. V. Mecking, N. S. Keenlyside, and R. J. Greatbatch. Multiple timescalesof stochastically forced north atlantic ocean variability: A model study.Ocean Dynamics, 65(9-10):1367–1381, 2015. → page 42[24] J. J. More´ and D. C. Sorensen. Computing a trust region step. SIAM Journalon Scientific and Statistical Computing, 4(3):553–572, 1983. → page 30[25] G. R. H. Paul Blanchard, Robert L. Devaney. Differential equations pp.96111. Cole Thomson Learning, 2006. → page 2[26] H. Peter D. Kernel estimation of a distribution function. Communications inStatistics-Theory and Methods, 14(3):605–620, 1985. → pages 29, 44[27] P. Petersen. Hysteresis in the bifurcation structure of a nonlinear opticaldevice with a time-dependent bifurcation parameter. In Electro-optic andPhotorefractive Materials, pages 371–376. Springer, 1987. → page 4[28] A. S. Pikovsky and J. Kurths. Coherence resonance in a noise-drivenexcitable system. Physical Review Letters, 78(5):775, 1997. → page 82[29] R. A. Serway, C. J. Moses, and C. A. Moyer. Modern physics pp. 212-213.Cengage Learning, 2004. → pages 32, 36, 46[30] J. Sieber and J. M. T. Thompson. Nonlinear softening as a predictiveprecursor to climate tipping. Phil. Trans. R. Soc. A, 370(1962):1205–1227,2012. → pages 18, 22, 23, 25, 31, 59, 61, 81, 94, 95, 97[31] S. Solomon, D. Qin, M. Manning, Z. Chen, M. Marquis, K. B. Averyt,M. Tignor, H. L. Miller, et al. Contribution of working group i to the fourth100assessment report of the intergovernmental panel on climate change ch.5,2007, 2007. → page 7[32] H. Stommel. Thermohaline convection with two stable regimes of flow.Tellus, 13(2):224–230, 1961. → page 5[33] S. Thorpe. The turbulent ocean, 439 pp. Cambridge University, 2005. →page 16[34] E. R. Widiasih. Dynamics of the budyko energy balance model. SIAMJournal on Applied Dynamical Systems, 12(4):2068–2092, 2013. → page 97[35] X. Wu and M. Mu. Impact of horizontal diffusion on the nonlinear stabilityof thermohaline circulation in a modified box model. Journal of PhysicalOceanography, 39(3):798–805, 2009. → pages 12, 16101Appendix ASupporting MaterialsA.1 Noise amplitude and drift rate chosen analysis fordifferent modelsIn this paper, we have investigated three different saddle node bifurcations in differ-ent models, and one may notice that we change the values of the noise amplitudesand drifting rates in different models. The criterion for choosing compatible noiseamplitude and drifting rate is to make sure these models can be reconstructed intoa same standard saddle node normal form.The standard saddle node normal form with additive noise isx˙ = µ+ x2+σ1dWt (A.1)And that is the same as Eq 3.1.In the 1D hysteresis model, saddle node is located on the upper branch:y˙ = (−υ+1)− (y−1)2+σ2dWt (A.2)Let υ1 =−υ+1 and y1 = y−1, then Eq A.2 is equivalent toy˙1 = υ1− y21+σ2dWt (A.3)102Since Eq (A.3) is exactly the standard saddle node form, we should use the samenoise amplitude and drifting rate in A.1 and A.2.To compare the Stommel model with the saddle node normal form, first we lin-earize the system. The system isdTdt= η1−T (1+ |Ψ|)+σ3dWtdΨdt= η1−η2−Ψ|Ψ|−T +(T −Ψ)η3+σ3dWtThe deterministic versions of Ψ and T are (consider T as a function of Ψ, and Ψas a function of T,η2):f (Ψ,η2) = Ψ˙= η1−η2−Ψ2−T (Ψ)+η3(T (Ψ)−Ψ)g(T,η2) = T˙ = η1−T (1+ |Ψ|)And when the system reaches the smooth bifurcation (Ψ∗, η∗2 ), f (Ψ∗,η∗2 ) = 0,g(T ∗,η∗2 ) = 0. From the cubic function root formula, we need η∗2 satisfies belowequation, which is the requirement for having two same real roots and one imagi-nary root.((1+η3)(η∗2 +η3−η1)6− (1+η3)327− η∗2 −η1η32)2 =−(η∗2 +η3−η13− (1+η3)29)3Now we expand the original equations around (Ψ∗, η∗2 ) and (T ∗, η∗2 ),Ψ˙=−(η2−η∗2 )+(−1+η3−12TΨΨ(Ψ∗))(Ψ−Ψ∗)2+σ3dWtT˙ = (−1−Ψ∗−T ∗ΨT (T ∗))(T −T ∗)+ 12(−2ΨT (T∗)−T ∗ΨT T (T ∗))(T −T ∗)2+σ3dWtlet β = η2−η∗2 , P = Ψ−Ψ∗, and re-scale P and β with A, where A = (−1+η3−12 TΨΨ(Ψ∗)), we getP˙ =−Aβ +P2+A ·σ3dWtIn this form, we can now compare it to the normal form A.1, and conclude that the103smooth bifurcation in the Stommel model is a saddle node. As for the values of theparameters, σ1 = A ·σ3, the drifting rate in A.1 is also A times larger than in theStommel model.104
Thesis/Dissertation
eng
Mathematics
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Noise induced tipping in ocean circulation based on the Stommel model
Text