Non-Smooth Dynamics in theStommel Model of the ThermohalineCirculationbyCody GriffithB.Sc., Metropolitan State University of Denver, 2016A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF SCIENCEinThe Faculty of Graduate and Postdoctoral Studies(Mathematics)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)August 2018c© Cody Griffith 2018The following individuals certify that they have read, and recommend to the Faculty ofGraduate and Postdoctoral Studies for acceptance, a thesis entitled:Non-Smooth Dynamics in the Stommel Model of the Thermohaline Circulationsubmitted by Cody Griffith in partial fulfillment of the requirements of the degree ofMaster of Science in Mathematics.Examining Committee:Rachel Kuske, Professor of MathematicsSupervisorBrian Wetton, Professor of MathematicsSupervisory Committee MemberiiAbstractWe analyzed the non-smooth dynamics for the Stommel model for thermohaline circulationwith additional mechanisms like slowly varying bifurcation parameters and high frequencyoscillatory forcing. Our goal was to find an analytic approximation to the tipping point andforced bifurcation induced by the new features in the model. We first analyze a simpler onecomponent model that has similar structure to the Stommel model and gradually build inmore complexity into the one component problem and study the effects on the tipping point.In this context, we compare the relative strengths of the non-smooth effects on the tippingpoint to the smooth effects which has been previously studied. With the one componentmodel understood, we then apply similar methods to the Stommel model and study theeffects on the non-smooth tipping points and forced bifurcations. With these results wehave the ability to fully describe the hysteresis found in the Stommel model.iiiLay SummaryWe study the behavior of the thermohaline circulation. This circulation is responsiblefor moving water from around the globe and thus a model was created to understandhow it functions. Work has previously been done on certain components of this model,but the analysis we provide is done around the less studied pieces known as the non-smooth dynamics. We ultimately provide a solution to the non-smooth dynamics and evenincorporate more mechanisms in the original model to account for a larger class of observablebehavior. This allows for a better understanding of the thermohaline circulation and couldhelp predict and prepare for sudden abrupt changes to the ocean currents that drive ourclimate.ivPrefaceThis thesis is original, unpublished, independent work by the author, C. Griffith.vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 One Component Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112.1 Static Bifurcations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112.2 Slowly Varying Bifurcation Parameter . . . . . . . . . . . . . . . . . . . . . 122.2.1 Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162.3 High Frequency Oscillatory Forcing . . . . . . . . . . . . . . . . . . . . . . 172.3.1 Case I: v0(t) ≤ −|A| . . . . . . . . . . . . . . . . . . . . . . . . . . . 202.3.2 Case II: |v0(t)|< |A| . . . . . . . . . . . . . . . . . . . . . . . . . . . 212.3.3 Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 242.4 Slowly Varying and Oscillatory Forcing . . . . . . . . . . . . . . . . . . . . 252.4.1 Case I: λ ≤ 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 272.4.2 Case II: λ > 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 292.4.3 Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 363 Two Component Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 393.1 Slowly Varying Bifurcation Parameter . . . . . . . . . . . . . . . . . . . . . 403.2 High Frequency Oscillatory Forcing . . . . . . . . . . . . . . . . . . . . . . 453.2.1 Case I: P0(t) ≤ −|A| . . . . . . . . . . . . . . . . . . . . . . . . . . . 483.2.2 Case II: |P0(t)|< |A| . . . . . . . . . . . . . . . . . . . . . . . . . . . 493.2.3 Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 533.3 Slow Variation with Oscillatory Forcing . . . . . . . . . . . . . . . . . . . . 553.3.1 Case I: λ ≤ 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 573.3.2 Case II: λ > 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 603.3.3 Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68viTable of Contents4 Summary and Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . 73Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79AppendicesA The Stommel Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81B One Component . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83C Two Component . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85viiList of Figures1.1 Vector field of a saddle-node bifurcation a = 0. . . . . . . . . . . . . . . . . 21.2 Bifurcation diagram of saddle-node bifurcation a = 0. . . . . . . . . . . . . 21.3 The saddle-node bifurcation. In (a) an example of tipping occurring aroundthe bifurcation for two sizes of slow variation, = {.01, .1}. In (b) a zoomin closer to the bifurcation. The dashed ( = .01) and dash-dotted ( = .1)black lines are the numerical solutions, we overlay the bifurcation diagramfor reference. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41.4 The Stommel Two Box Model: Differing volume boxes with a temperatureand salinity, Ti and Si. The boxes are connected by an overflow and capillarytube that has a circulation rate V . There is also a surface temperature andsalinity for each box, Tis and Sis. We assume that there is some stirring togive a well mixed structure. . . . . . . . . . . . . . . . . . . . . . . . . . . . 61.5 The equilibria of the non-dimensionalized system (1.6). Parameters valuesare η1 = 4 and η3 = .375. The above plots are two-dimensional projectionsof the full 3-dimensional system in (η2,V ,T ). We see non-smooth behaviorhappening in both plots when V = 0. The red line indicates a stable branchwhere the dashed dotted line is for an unstable branch. . . . . . . . . . . . . 71.6 The choice in η3 dictates the orientation of the problem, in each plot we havefixed η1 = 4. The case for η3 = 1 is special due to the two bifurcationsoverlapping and the unstable equilibrium vanishing. . . . . . . . . . . . . . 82.1 The one component bifurcation diagram with the upper and lower equilib-rium branches as well as the unstable middle branch. The non-smooth bi-furcation occurs at (0,0) denoted by the circle and the smooth bifurcationoccurs at (1,1) by the box. . . . . . . . . . . . . . . . . . . . . . . . . . . . 122.2 In (a) the numerical solutions (black dashed and dash-dotted lines) to (2.1)are given with A = 0 and = {.01, .04} respectively. The bifurcation plotis overlayed for convenience. In (b) a zoom in of what happens near thenon-smooth bifurcation. The solid vertical lines (black) are tipping pointswhere we use the tipping criterion x > .5 on the numerical solution. Thedashed and dash-dotted vertical lines (blue) are the tipping estimates. In (c)a range of and their corresponding tipping (red stars) are compared to ourestimate (solid black line) from (2.14). . . . . . . . . . . . . . . . . . . . . . 152.3 The parameter ranges for each case are shown here with A = 2. For reference,the original bifurcation diagram is overlayed. . . . . . . . . . . . . . . . . . 20viiiList of Figures2.4 The non-smooth function |y0(T )|= |v0−A cos(T )| that we integrate is shownas a solid red line. We also show an example of v0 as a horizontal blue dottedline. Here the value of |v0|≤ |A|, which causes kinks to appear at the rootsof |y0|: T1 and T2 respectively. These are the vertical black dashed dottedlines. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 212.5 In (a) the numerical time series solutions to (2.1) are given from bottom to topwith µ = {.8, .33, .15} in case I, case II and µ < µosc in (2.41) respectivelywith A = 2, Ω = 10 and = 0. In (b) we show the time series on thebifurcation diagram. In (c), a zoom in closer to the non-smooth bifurcationof (b), where the dotted vertical lines dictate the region between case I andcase II (green) as well as the bifurcation estimate (blue) respectively. In (d)a range of Ω−1 and the corresponding numerical bifurcations (red stars) arecompared to our estimate of the bifurcations (black solid line). We considerthe bifurcation criterion to be when the numerical solution has passed x > .5. 232.6 On the left, one can see the bifurcation diagram for the canonical system(2.1) with the numerical solution (black dotted line). On the right, a zoomin around the non-smooth bifurcation. The dotted vertical line is the tippingpoint µmixed (2.66) (blue). The vertical line (black) is the when the numericalsolution has passed the tipping criterion x > .5. The parameter values are = .05, λ = .8 and A = 4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 322.7 On the left, one can see the bifurcation diagram for the canonical system(2.1) with the numerical solution (black dotted line). On the right, a zoom inabout the non-smooth bifurcation. The dotted vertical lines are the tippingpoint µmixed (2.66) (blue) and slowly varying tipping µslow (2.14) (green).The vertical line (black) is the when the numerical solution has passed thetipping criterion x > .5. The parameter values are = .05, λ = 1.05 andA = 4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 332.8 On the left, one can see the bifurcation diagram for the canonical system(2.1) with the numerical solution(black dotted line). On the right, a zoom inaround the non-smooth bifurcation. The dotted vertical lines are the tippingpoint µmixed (2.66) (blue) and slowly varying tipping µslow (2.14) (green).The vertical line (black) is the when the numerical solution has passed thetipping criterion x > .5. The parameter values are = .05, λ = 1.6 and A = 4. 342.9 An example of numerical tipping (red stars) as the numerical solution to (2.1)passes the tipping criterion x = .5 for the last time. The parameter valuesare = .01 and A = 4. The lines are the case I tipping estimate (2.66) (blacksolid line) and the case II tipping estimate (2.14) (blue dotted line). . . . . 352.10 The numerical tipping (red stars) follows the appropriate case depending onλ. The case I tipping estimate (black solid line) and the case II tippingestimate (blue dotted line) are shown. . . . . . . . . . . . . . . . . . . . . . 35ixList of Figures3.1 In (a) the numerical solutions (black dotted and dash-dotted lines) to (3.1)are given with η1 = 4, η3 = .375, and = {.01, .04} respectively. In (b) azoom is given closer to the non-smooth bifurcation region. The blue verticallines are the predictions (3.17) against the black solid vertical lines which arethe tipping points with the tipping criterion V > Vsmooth on the numericalsolution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 433.2 In (a) we have the numerical solutions (black dotted and dash-dotted) overthe standard equilibrium plot for V vs. T . In (b) a zoom in closer to thebifurcation area. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 443.3 The numerical tipping (red stars) vs. the estimate η2slow (black line) withη1 = 4 and η3 = .375. The top blue line is the tipping of the second eigen-value. The tipping criterion on the numerical solution is V > Vsmooth. . . . 453.4 Here we have the parameter ranges for case I and case II shown as the rightmost green vertical line and the bifurcation value as the left blue vertical linerespectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 483.5 In (a) the numerical time series solutions to (3.1) is given with parameters ineach qualitatively different case of η2 = {2.3, 1.8, 1.51} with η1 = 4, η3 = .375,A = B = 2 and Ω = 10. In (b) these same solutions are shown on thebifurcation diagram. In (c) a zoom in closer to the non-smooth bifurcationregion where the blue vertical line is the estimated bifurcation (3.41). . . . 513.6 In (a) we have the same numerical time series solutions for the qualitativelydifferent cases η2 = {2.3, 1.8, 1.51}. In (b) we plot these solutions over thestandard equilibrium plot for V vs. T . In (c) a zoom closer to the bifurcationarea is provided. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 523.7 The numerical tipping (red stars) vs. the estimate (black line). The modelparameters are η1 = 4,η3 = .375 and A = B = 2. The bifurcation criterionfor the numerical solution is V > Vsmooth. . . . . . . . . . . . . . . . . . . . 533.8 The model values are λ = .8, = .01 with A = B = 2. In (a) the numericalsolution (black dotted line) to (3.1) is given with η1 = 4, η3 = .375. In(b) a zoom in closer to the non-smooth bifurcation region where the bluedotted vertical line is the tipping point (3.70) and the black vertical line arethe tipping points with the tipping criterion V > Vsmooth on the numericalsolution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 643.9 The model values are λ = .8, = .01 with A = B = 2. In (a) we have thenumerical solution (black dotted) over the static equilibrium plot for V vs.T . In (b) a zoom of the bifurcation area. . . . . . . . . . . . . . . . . . . . . 643.10 The model values are λ = 1.05, = .01 with A = B = 2. In (a) the numericalsolution (black dotted line) to (3.1) is given with η1 = 4 and η3 = .375. In(b) a zoom in closer to the non-smooth bifurcation region where the bluedotted vertical line is the mixed tipping point (3.70), green dotted verticleline is the slow tipping point (3.17) and the black vertical line are the tippingpoints with the tipping criterion V > Vsmooth on the numerical solution. . . 653.11 The model values are λ = 1.05, = .01 with A = B = 2. In (a) we have thenumerical solution (black dotted) over the static equilibrium plot for V vs.T . In (b) a zoom of the bifurcation area is given. . . . . . . . . . . . . . . . 65xList of Figures3.12 The model values are λ = 2, = .01 with A = B = 2. In (a) the numericalsolution (black dotted line) to (3.1) is given with η1 = 4 and η3 = .375. In (b)a zoom in closer to the non-smooth bifurcation region where the blue dottedvertical line is the mixed tipping point (3.70), the green dotted vertical lineis the slow tipping point (3.17) and the black vertical line are the tippingpoints with the tipping criterion V > Vsmooth on the numerical solution. . . 663.13 The model values are λ = 2, = .01 with A = B = 2. In (a) we have thenumerical solution (black dotted) over the static equilibrium plot for V vs.T . In (b) a zoom of the bifurcation area is provided. . . . . . . . . . . . . . 663.14 An example of numerical tipping (red stars) as the numerical solution to (3.1)passes the tipping criterion V > Vsmooth. The parameter values are = .01and A = B = 3. The lines are the case I tipping estimate (black solid line)and the case II tipping estimate (blue dotted line). . . . . . . . . . . . . . . 673.15 The numerical tipping (red stars) follows the appropriate case depending onλ for = 0.01. The case I tipping estimate η2mixed (black solid line) andslowly varying tipping estimate η2slow (blue dotted line) are shown. . . . . . 684.1 Comparison between the slow tipping points across . The blue solid lineis the non-smooth tipping points where the red dotted line is the smoothtipping points. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 744.2 Comparison between the oscillatory bifurcation across Ω−1 for A = 2. Theblue solid line is the non-smooth case where the red dotted line is the smoothcase. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 754.3 Comparison between the mixed tipping in (a) with a fixed Ω−1 = .1 and (b)with a fixed = .1. The blue solid line is the non-smooth case where the reddotted line is the smooth case. . . . . . . . . . . . . . . . . . . . . . . . . . 764.4 Low Frequency: Model parameters are = .01, A = B = 1 and Ω = 3. . . . 774.5 Large Amplitude: Model parameters are = .01, A = B = 300 and Ω = 1000. 774.6 Stochastic: In (a) many realizations of the numerical solution for V in (4.1)are given with model parameters η1 = 4, η3 = .375, = .01 and A = B = .7.In (b) a zoom in closer to the non-smooth bifurcation region can be seen. In(c) we have the realizations over the standard equilibrium plot for V vs. T .In (d) a zoom of the bifurcation area is shown. . . . . . . . . . . . . . . . . 78xiAcknowledgmentsI would like to thank my supervisor Rachel Kuske for the project idea and guidance as well asmy committee member Brian Wetton for helping me professionally through graduate school.Special thanks are given to my peers Tim Jaschek and Matthias Klo¨ckner for their aide inwriting and support throughout the entire degree. Couldn’t have done it without you both!xiiChapter 1IntroductionDynamical systems is the study of the possible states an observable solution may experienceand is important in most engineering, biological or even chemical systems to name a few.This approach allows conditions to be given for when a solution can be found or whenthere is stable behavior. Often we find that parameters inherent in the model play hugeroles in the dynamical behavior and they can be the difference between a system having anequilibrium or not. When we find a parameter that has this effect, we call it a bifurcationparameter since there is some value that changes the qualitative behavior of the system.For example, the Hodgkin-Huxley model for neurons contains a parameter for injectedcurrent I which turns out to be a bifurcation parameter with a Hopf bifurcation. Thisbifurcation is responsible for the actual firing of a neuron in the brain. In epidemiologicalmodeling, the basic SIR model with an additional transition function between the infectedand recovered population causes the reproduction number R0 to become a bifurcation pa-rameter with a backward bifurcation. This causes a temporary equilibrium to form in theinfected population that usually would never see an equilibrium. Even in activation poten-tials of neural networks, using a hyperbolic tangent function causes a bifurcation to occurin the synaptic feedback parameter w which results in a pitchfork bifurcation. This causeswildly different equilibria for learned parameters in a machine learning setting. The canon-ical example is the saddle-node bifurcation and was the first to be found within a complexsystem studied from a dynamical perspective. The saddle-node bifurcation has the locallytopological equivalent formx˙ = a− x2, (1.1)where by locally topological equivalence we mean the behavior near the bifurcation may berepresented in this form. This property is critical to reducing most complex problems andmodels to simpler local problems that can be studied individually.1Chapter 1. IntroductionFigure 1.1: Vector field of a saddle-node bifurcation a = 0.Figure 1.2: Bifurcation diagram of saddle-node bifurcation a = 0.In figure 1.1 we show the vector field of the system that contains a saddle-node bifurca-tion. The equilibria of this system is x = ±√a for any a ≥ 0 where stable equilibrium pointsare marked with red filled points and unstable with unfilled points. Notice that at a = 0we are no longer able to find a stable equilibrium and when a < 0 there are no equilibria atall. Thus a = 0 is a simple example of a bifurcation where two equilibria annihilate. This2Tipping in a Slowly Varying Systembehavior is why the bifurcation is often referred to as a fold bifurcation, although we referto this as a saddle-node bifurcation in this thesis. In figure 1.2 we plot the same systemagainst the parameter a, which we call the bifurcation diagram. Here we see the regionwith two equilibrium a > 0, the bifurcation a = 0 and the region of no stability a < 0. Formore on the saddle-node bifurcation see [12].There are many types of bifurcations that appear in different systems that each havetheir own key properties. Studying these properties leads to a deeper understanding of thesystem on both a global and a local scale. Work has been done on systems that have smoothbifurcations due to how commonly these appear, but non-smooth dynamics still are presentin the physical world.Non-smooth bifurcations are a topic that arise in special systems and for how frequentthey appear, they have not been studied nearly as much as their smooth counter parts.This thesis discusses the role of the non-smooth saddle-node bifurcation in a simplified onecomponent system in chapter 2 as well as in the classic Stommel model for thermohalinecirculation dynamics in chapter 3. Many interesting ocean and weather mechanisms maybe incorporated into the Stommel model to provide more realistic predictions for weatherpatterns. We choose to study slowly varying bifurcation parameters and their effect on thestability of a system while contrasting this with non-autonomous oscillatory forcing. Theinteraction of these features causes complex dynamics around the standard bifurcations andcan lead to an advanced bifurcation or delayed tipping. For the one component system, adetailed analysis of these features is done on the smooth bifurcation in [24].Tipping in a Slowly Varying SystemA system with a parameter known to cause a bifurcation will no longer admit a bifurcationin the standard sense when the parameter slowly varies. Instead, these conditions give riseto a smooth but rapid change in the system’s equilibria. The point in which this behavioroccurs is then called a tipping point.More formally, a tipping point is a point that causes an abrupt smooth transition indynamical behavior as the system moves into a qualitatively different state. This is usuallycaused by some exterior control system that pushes change towards a different state oncea critical point has been passed, for example with biological systems seen in [2]. Theseare known to be caused by changes in one or more parameters in the system. An analysisthat lays the theoretical backing of slowly varying parameters with algebraic bifurcationsis found in [7].Tipping points have been discovered to occur in a wide variety of systems and have be-come a big staple in the study of areas like catastrophe theory and dynamical systems. Theyaid in predicting the future of a system and even could be a warning for irreversible changelike in the case of the Stommel model. A tipping point thus shares similar characteristicsto a bifurcation and typically occurs close to the static bifurcation location.In this thesis we use the results from [24] where the systemx˙ =Da+ k0 + k1x+ k2x2,a˙ =− , (1.2)3Tipping in a Slowly Varying Systemwhere 1 was studied. This model is a slowly varying quadratic differential equationcontaining a smooth saddle-node bifurcation and appears in many physical models, forexample [6]. A key result from [24] is that the solution and tipping point for (1.2) have theformx ∼ 1|k2|(k12+(D|k2|)1/3) Ai′ ([D|k2|]−2/3 (k214 + k0|k2|+D|k2|a))Ai([D|k2|]−2/3(k214 + k0|k2|+D|k2|a)) (1.3)atip = (D|k2|)−1/3aAiry − asDfor as = k0 +k214|k2| , (1.4)with Ai(·) being the Airy function and aAiry = 2/3 · (−2.33810 . . .) corresponding to thefirst zero of the Airy function. The singularity found in (1.4) is a recurring tool for the workpresented in this thesis, even though we deal with a version of (1.2) that has a non-smoothbifurcation.(a) (b)Figure 1.3: The saddle-node bifurcation. In (a) an example of tipping occurring aroundthe bifurcation for two sizes of slow variation, = {.01, .1}. In (b) a zoom in closer to thebifurcation. The dashed ( = .01) and dash-dotted ( = .1) black lines are the numericalsolutions, we overlay the bifurcation diagram for reference.In figure 1.3 we show a numerical solution to the simple saddle-node system with tipping(1.2). Here we have D = 1, k0 = k1 = 0 and k2 = −1 which is the model from (1.1). Thesolution follows closely to the stable branch even after the bifurcation for the static model,which is an example of this delayed behavior. From here on we refer to the numerical tippingpoint to be when the numerical solution has passed a threshold away from the equilibriumsuch that it is reasonable to say the solution is transitioning to a new state. We call thisthreshold the tipping criterion and it is specified whenever we compare our estimates to thenumerical solution.The task of finding where tipping points occur depends on the situation, but in generalthe approach is to search for when a solution to a model fails or becomes large. Examples4The Stommel Modelcould be when the solution fails to be real or when an exponential term grows large, bothof which are seen throughout this thesis.The Stommel ModelGlobal circulation models have primarily focused on three different categories:• Atmospheric components - the effect greenhouse gases have on the atmosphere,• Oceanic components - the effect of tides and interaction of temperature with salinityin the oceans,• Sea ice and land surface components.These categories all contribute significantly to the overall prediction of weather and climatefor the planet, which has importance to just about every industry and economy. Failure toadhere to and prepare for sudden changes in the climate has led to drastic situations likesevere droughts or ocean acidification. Atmospheric models have been vastly studied butfar less work has been done on the contribution from the ocean and the dynamics that drivethe tides and currents.A key feature of oceanic models is when patterns form around regions of bi-stability oftemperature and salinity. An example of this is the thermohaline circulation (THC) whichhas abrupt qualitative changes at certain points, see [1, 14, 17, 18]. Just earlier this yearevidence was found of weakening occurring around these abrupt changes in a system ofocean patterns known as the Atlantic meridional overturning circulation (AMOC) [4]. Thisis the first evidence of ocean dynamics responding to temperature change on the surfaceand can help further predict the future of the system. It is imperative that appropriateaction is taken to prepare for the future of these type of systems as they are outside ourrealm of control.To study these phenomena we create parametric models to replicate the dynamics weobserve. Initially, Henry Stommel proposed the two box model in 1961 to understand thephysics of the THC, shown in figure 1.4. In [23], it is suggested that there are actuallytwo different stability regimes which even overlap in the system that is proposed and con-cluded that oceanic dynamics behave very similarly about these equilibria. These type ofsystems have since been a heavily studied area for both climatology due to the wide rangingapplications and dynamical systems for its generalization into bi-stability.5The Stommel ModelFigure 1.4: The Stommel Two Box Model: Differing volume boxes with a temperature andsalinity, Ti and Si. The boxes are connected by an overflow and capillary tube that has acirculation rate V . There is also a surface temperature and salinity for each box, Tis andSis. We assume that there is some stirring to give a well mixed structure.With emphasis on mathematics, the focus of this thesis is on developing an effectiveapproach to models with bi-stability and additional mechanisms. Thus the physical quan-tities are brushed aside in favor of their non-dimensional alternatives, see Appendix A forthe derivation. The non-dimensionalized Stommel model is represented with the systemT˙ = η1 − T (1 + |T − S|),S˙ = η2 − S(η3 + |T − S|).(1.5)The variables T and S are the temperature and salinity respectively where the non-smoothnessis seen directly from the |T − S| term. The parameters η1, η2, and η3 are all dimensionlessquantities that each have physical interpretation to the relaxation times and volumes of thebox. Here η1 is thought of as the thermal variation, η2 as the saline variation otherwiseknown as the freshwater flux, and η3 as the ratio of relaxation times of temperature andsalinity. It also is a physical restriction for both η1 and η3 to be positive quantities thattake any value. The parameter η3 has the additional property to determine the orientationof the equilibria. We denote a standard orientation to be when η3 < 1, reverse orientationfor η3 > 1, and η3 = 1 a special case. The different orientations are shown in figure 1.6.Recall that η3 is the ratio of relaxation rates and when η3 = 1 the relaxation rates for boththe thermal and salinity variables are the same. Under these conditions we lose bi-stabilityand instead see a single stable equilibrium.6The Stommel ModelThe parameter η2 is the most interesting as different values cause major qualitativeand quantitative changes in the dynamics of the system. Bifurcations have been discoveredat two different points in the system, each being called either a smooth or a non-smoothsaddle-node bifurcation. In the Stommel model, it is convenient to view the system in termsof the circulation rate V = T − S, see Appendix A for the derivation. This leads to thesystemT˙ = η1 − T (1 + |V |),V˙ = (η1 − η2)− V |V |−T + η3(T − V ).(1.6)(a) V vs. T (b) η2 vs. VFigure 1.5: The equilibria of the non-dimensionalized system (1.6). Parameters valuesare η1 = 4 and η3 = .375. The above plots are two-dimensional projections of the full3-dimensional system in (η2,V ,T ). We see non-smooth behavior happening in both plotswhen V = 0. The red line indicates a stable branch where the dashed dotted line is for anunstable branch.As shown in figure 1.5, the equilibrium curves reveal much about the dynamics. In (a)the graph of the equilibria for V vs. T shows non-smooth behavior occurring at V = 0 andin (b) the two types of bifurcation appear clearly in the graph of equilibria for η2 vs. V . Inthis plot, both the upper and lower branches of the equilibrium are stable with the middlebranch being unstable. The stable branches relate to which variable is dominant. For thelower branch, we call this the saline branch, and the upper branch the thermal branch. Thelocation of the non-smooth bifurcation is found analytically, (η2ns, Vns, Tns) = (η1η3, 0, η1),and the smooth bifurcation, (η2smooth, Vsmooth, Tsmooth), is the only real solution to a cubicpolynomial. The smoothness of each bifurcation is apparent and arises from the absolutevalue term in the defining dynamics of (1.6), which is non-smooth only at V = 0.7The Stommel Model(a) η3 = .375 (b) η3 = 1(c) η3 = 1.875Figure 1.6: The choice in η3 dictates the orientation of the problem, in each plot we havefixed η1 = 4. The case for η3 = 1 is special due to the two bifurcations overlapping and theunstable equilibrium vanishing.Much is known about the Stommel model in the case where η2 is fixed but realisticallythis is not the case. In [17], this parameter is described as the influx of freshwater into theAtlantic and the changing nature of η2 is justified by a positive feedback loop for salinitythat drives the THC to move high-salinity water towards deep oceans. This loop causes theabrupt smooth bifurcation but then afterwards, a salinity deficit causes the parameter todecrease back towards the non-smooth bifurcation.This type of behavior is known as hysteresis, where there is some bi-stability region thatthe solution cycles through and observes both states of the equilibria. A similar analysis tothe Stommel model’s hysteresis can be found in [19]. The phenomena of hysteresis appearsin many physical systems, for example [11, 8, 10]. The smooth component of the hysteresiscurve has been studied in a reduced one component model, see [24]. In this thesis wecomplement these results with an analysis of the non-smooth component.8Numerical MethodsNumerical MethodsTo obtain numerical solutions to the ordinary differential equations studied in this thesis, wechoose to use both 2nd order and 4th order Runge-Kutta methods. The 2nd order methodare used for the one component model and the 4th order method for the two componentmodel. The choice in these methods comes from using the simplest scheme since numericalsensitivity is not present in our problem. We use the numerical solutions to compare ourapproximations to the observed state transition.9Chapter 2One Component ModelWe consider a simpler system to give insight into the more complex two component Stommelmodel. Here we use a toy system to build the analysis on and the spatial variables may nothave a physical interpretation. This system is the following one component model in termsof the variables x and µx˙ = −µ+ 2|x| − x|x|+A sin(Ωt),µ˙ =− , (2.1)x(0) = x0, µ(0) = µ0,where the fixed parameters are the slow variation rate of 1, the amplitude of oscil-lation A and the frequency of oscillation Ω. We also assume the initial conditions to bex0 = 1−√1 + µ0 and µ0 > µns which focuses our calculations on the lower equilibriumbranch where x < 0 and study nearby behavior. The value µns refers to the non-smoothbifurcation which is discussed below in section 2.1.The system (2.1) is generalized from a basic model that contains both a smooth andnon-smooth saddle-node bifurcation. This structure is similar to the Stommel model andhence a good model to test features like slow variation or oscillatory forcing. The slowvariation is clear, but we use oscillatory forcing here in preparation for the two componentmodel of the next chapter. In each case, emphasis is put on the non-smooth component ofthe model to study the non-smooth bifurcation and the role it plays in the hysteresis curvewe anticipate in the Stommel model.2.1 Static BifurcationsThe foundation to our understanding comes from the simplest structure lying within thecanonical system (2.1) which is the bifurcation structure. This means finding the generalform for the equilibria in (2.1) with A = 0 and = 0, which is our basic model with a staticµ and no forcing. As we have a fixed parameter value, we search for a point or set of pointsthat the solution relaxes to as t→∞. We call these points the equilibrium points and theyare either stable or unstable. Since we are considering all possible µ, we want all of theequilibrium points for each µ and thus we call these the equilibrium branches.To find all equilibrium branches, we search for when the solution has come to a rest,which is equivalent to setting the derivative of x to zero. Thus we set (2.1) to zero with0 = −µ+ 2|x|−x|x|. (2.2)Solving (2.2) results in 3 solutions where the stability of each is characterized by smallperturbations to the equilibrium either linearly growing or decaying. We denote the stable102.2. Slowly Varying Bifurcation Parameterequilibria as xl and xu for the lower and upper branches respectively, and a single unstablemiddle branch, xm. These are given byxl = 1−√1 + µ, xu = 1 +√1− µ, xm = 1−√1− µ.We note that xl is valid for µ ≥ 0 and both xu and xm for µ ≤ 1. Thus this systemhas a stable equilibrium for each value of the parameter and has a region of bi-stability for0 ≤ µ ≤ 1. The boundaries of this region are (µns, xns) = (0, 0) and (µsmooth, xsmooth) =(1, 1) which are the non-smooth and smooth saddle-node bifurcations respectively. Bothare saddle-node due to pairs of equilibria annihilating at these locations which is shown infigure 2.1.Figure 2.1: The one component bifurcation diagram with the upper and lower equilibriumbranches as well as the unstable middle branch. The non-smooth bifurcation occurs at (0,0)denoted by the circle and the smooth bifurcation occurs at (1,1) by the box.2.2 Slowly Varying Bifurcation ParameterTo develop a method for the slowly varying Stommel model, we consider (2.1) with 1and A = 0. Under these conditions, µ(t) is a function of time and thus a bifurcationno longer occurs. Instead, it is expected that a tipping point occurs nearby the staticbifurcation points as long as is small. Also, due to µ(t) being a function of time, we willfind equilibria that are also functions in time, which we call pseudo-equilibira. The smoothcase is well understood, see [24], so we consider the behavior of the non-smooth bifurcationwith x < 0. From [7] as well as the smooth model [24], it is common practice to rescaletime in a model with slow variation to put the dynamics on the same order and allow foralgebraic solutions to be found. Here the parameter µ(t) is slowly varying in time so it112.2. Slowly Varying Bifurcation Parametermakes sense to rescale using this as our slow time, τ = t. Applying both x < 0 and thisslow time approach to the system (2.1) then givesxτ =− µ(τ)− 2x+ x2,µτ =− 1.(2.3)A standard approach to extracting information out of complicated models is to findreduced equations by separating the behavior at each order of the slow time. This approachis known as using an asymptotic expansion and further details can be found in Murray’sAsymptotic Analysis [15]. With being the small quantity that dictates our slow time, wechoose to use an asymptotic expansion of x withx(τ) ∼ x0(τ) + x1(τ) + 2x2(τ) +O(3). (2.4)This approach captures the slowly varying behavior of the solution in terms of this smallquantity and aims to relate the slow variation to the solution. We substitute the expansion(2.4) into the scaled system (2.3) to getx0τ + 2x1τ + . . . = −µ(τ)− 2x0 + x20 + (−2x1 + 2x1x0) + 2(−2x1 + 2x2x0 + x21) + . . .Once we separate the equations at each order of , we find the following system of equationsO(1) : 0 = −µ(τ)− 2x0 + x20, (2.5)O() : 0 = −x0τ − 2x1 + 2x1x0, (2.6)O(2) : 0 = −x1τ − 2x2 + 2x2x0 + x21. (2.7)Each of the equations (2.5)-(2.7) gives the respective order’s pseudo-equilibrium. Thus wesolve each equation progressively to find the terms of our asymptotic expansion (2.4) asx(t) ∼ 1−√1 + µ(t) +4(1 + µ(t))− 3232(1 + µ(t))5/2+O(3). (2.8)We call (2.8) the outer solution as it approximates the solution well for values of x(t) awayfrom the bifurcation value µns. Since the dynamics of the system (2.1) change at x = 0 dueto the non-smooth bifurcation of the underlying static system, this solution is valid only forx < 0 and µ > 0.It is a key assumption of an asymptotic expansion that the terms are clearly separatedby order of . We search for a scaling of µ and x for which (2.8) is no longer valid under thisassumption of order separation. This assumption fails when x0 ∼ x1 which occurs herefor µ ∼ O(). To confirm, we conduct a simple scale analysis to determine the appropriatescaling for the local analysis about x = 0. Hence we consider the general scalesx = αy, µ = βm,122.2. Slowly Varying Bifurcation Parameterwith α > 0 and β > 0 for an inner scaling. We apply these local variables in (2.1) to getthe systemαy˙ =− βm+ α2|y|−2αy|y|,βm˙ =− . (2.9)We balance the leading order terms αy˙ with βm to find α = β. Here the equation for mcalls for β = 1, thus we have the scaling for the local analysisx = y, µ = m. (2.10)We have found that the scalings in (2.10) apply to all x and thus we consider the regionof x > 0. Substituting the local variables (2.10) into the original model (2.1) we find thefollowing inner system for the region of x > 0y˙ =−m(t) + 2y − y2,m˙ =− 1. (2.11)We recall that we are searching for a link between y and m, and from [7] we use that it isthen convenient to change the differentiation on y to be with respect to the slowly varyingparameter m. This incorporates the behavior of m(t) directly into the equation we solveand gives us a direct method for finding the tipping point. Then the leading order equationisym = m− 2y. (2.12)The leading order solution to (2.12) is found explicitly as followsy(m) = Ce−2m +m2− 14+O().With the inner solution found in terms of the parameter m, we write this in terms of theoriginal variables withx(t) ∼ Ce−2µ(t)/ + µ(t)2+O(). (2.13)We call the solution (2.13) the inner solution as it approximated the solution well near thebifurcation value µns. Since the inner solution behaves exponentially, the tipping point,µslow, occurs when the exponential term begins to grow rapidly. Here we consider tippingto occur when the solution becomes O(1/). Then we find the tipping point µslow to takethe formµslow =12 log(). (2.14)Thus we have the tipping point for the slowly varying model. Notice that for smallvalues of , µslow < µns and this is consistent with considering the inner equation (2.11) forthe region x > 0 as we found in the analysis. Hence we find that a slowly varying bifurcationparameter causes a delay in the rapid transition to the upper branch and we expect the132.2. Slowly Varying Bifurcation Parametersolution to remain near the lower branch for longer than in the static problem. In terms ofhysteresis, then slow variation allows for a longer period before the states switch from thelower to the upper branch.In figure 2.2 (a,b), two examples of this tipping is shown for different sizes of alongwith the standard bifurcation diagram where (c) demonstrates the tipping approximationacross a range of . The concavities match as well as agreement in the estimation of thetipping point as goes to 0.(a) (b)(c)Figure 2.2: In (a) the numerical solutions (black dashed and dash-dotted lines) to (2.1)are given with A = 0 and = {.01, .04} respectively. The bifurcation plot is overlayedfor convenience. In (b) a zoom in of what happens near the non-smooth bifurcation. Thesolid vertical lines (black) are tipping points where we use the tipping criterion x > .5 onthe numerical solution. The dashed and dash-dotted vertical lines (blue) are the tippingestimates. In (c) a range of and their corresponding tipping (red stars) are compared toour estimate (solid black line) from (2.14).142.2. Slowly Varying Bifurcation Parameter2.2.1 StabilityFrom the static model we know our outer solution (2.8) to be stable, but to verify thatthe inner solution (2.13) is stable we use a simple linear stability analysis on the innersystem. Typically to do this, an analysis would be performed about an equilibrium to seeif perturbations would grow or decay. Although, in this model there is a parameter thatis allowed to vary and hence we must be careful to note the analysis is about the pseudo-equilibrium instead. In the first region of interest, m(t) ≥ 0, the following inner equationand pseudo-equilibrium, z0(t), hold below the axisy˙ = −m(t)− 2y = f(t, y), z0(t) = −m(t)2. (2.15)We then consider simple perturbations of the pseudo-equilibrium, u, in (2.15) withy(t) = z0(t) + u(t), ‖u(t)‖ 1.Normally, a Taylor expansion would result in expressing the perturbations with their ownequation that we could use to determine stability. Since z0(t) is not fixed, we must considerits contribution to the derivative in this region of the parameter space with m(t) ≥ 0. Thuswe findy˙ = z˙0 + u˙,z˙0 ={−12m˙ = 12 m(t) > 0,0 m(t) = 0.(2.16)Now we apply the standard Taylor expansion to see the behavior of these perturbations andwith the contributions in (2.16), the inner equation (2.15) becomesy˙ =f(t, z0) + fy(t, z0)(y − z0) = fy(t, z0)u,u˙ ={−12 − 2u, m(t) > 0,−2u, m(t) = 0.(2.17)If this were the static parameter problem, we would always have the second case in(2.17), which is always stable due to the sign. Since we allow for a varying parameter,we learn that the solution is attracted to just below the pseudo-equilibrium z0(t). As thissystem always experiences the critical point m = 0 due to the continuous decrease in m(t),the slowly varying parameter eventually acts like the static parameter in section 2.1. Hencewe have that for x < 0, the pseudo-equilibrium is hyperbolic and asymptotically stable.Here there is a critical point at (µns, xns) = (0, 0) which corresponds to a non-hyperbolicequilibrium point. Generally, non-hyperbolic behavior signals equilibrium structures tochange. Here, this signals a transition in behavior for x > 0 and helps identify that thetipping occurs in this region.For the second region of interest, m(t) < 0, we found a solution that had the followinginner equation which has the pseudo-equilibrium above the axis withy˙ = −m(t) + 2y, z0(t) = m(t)2. (2.18)152.3. High Frequency Oscillatory ForcingIn (2.18) we find a contradiction, here m(t) < 0 yet the solution of this region is above theaxis x = 0. Thus we may conclude that this inner equation has no equilibrium in this regionand further verifies that the critical point (µns, xns) was non-hyperbolic and tipping occursfor m(t) < 0.2.3 High Frequency Oscillatory ForcingTo understand the oscillatory forcing in the Stommel model, consider the canonical system(2.1) with A ∼ O(1), Ω 1 and = 0, which gives high frequency oscillatory forcing in thesystem. Under these conditions, we have a static parameter and for each parameter valuethere is oscillatory forcing with solutions characterized by oscillations about a fixed point.Thus we should expect to find a bifurcation influenced by oscillations occurring under theseconditions. Here we develop a method to find oscillatory solutions to determine what theeffect of oscillatory forcing has on the bifurcation of (2.1). In section 2.2, we focused only onthe slowly varying dynamics but here we have both a slow time scale t and a fast time scaleT = Ωt. This naturally suggests a multiple scales approach where we search for a solutionthat is dependent on both of these scales, x(t) = x(t, T ). This method is commonly usedin problems that have behavior observable on multiple scales, and we use it here to finda way to accurately analyze each scale and effectively combine their behavior into a singleunifying solution. Further discussion on this method can be found in [20].Recall that our focus is on the non-smooth behavior and hence we restrict the solutionto follow along the lower stable equilibrium branch where x < 0. Using this multiple scalesapproach, our canonical system (2.1) has the following formxT + Ω−1xt = Ω−1(−µ− 2x+ x2 +A sin(T )) . (2.19)Note: We choose to use the subscript notation for partial derivatives, ∂x∂T = xT . In (2.19),the small quantity Ω−1 appears which suggests an asymptotic expansion in powers of thisquantityx(t, T ) ∼ x0(t, T ) + Ω−1x1(t, T ) + Ω−2x2(t, T ) +O(Ω−3). (2.20)Substituting (2.20) into (2.19), we findx0T + Ω−1x0t + Ω−1x1T + . . . = Ω−1(−µ− 2x0 + x20 +A sin(T )) + Ω−2(−2x1 + 2x1x0) + . . .Here we separate by each order of Ω to find the following system of reduced equationsO(1) : x0T = 0, (2.21)O(Ω−1) : x1T + x0t = −µ− 2x0 + x20 +A sin(T ), (2.22)O(Ω−2) : x2T + x1t = −2x1 + 2x0x1. (2.23)With an equation at each order, we must be able to solve each equation to proceed tothe next but we must also further restrict our solution from having resonant or linearly162.3. High Frequency Oscillatory Forcinggrowing terms to prevent any multiplicity or exponential growth. This assures that theterms in the asymptotic expansion are compatible with one another and we are able to finda robust solution. A common method to guarantee compatible solutions with sublineargrowth at each order is the Fredholm alternative. This provides a solvability condition foreach equation of the form xiT = Ri(t, T ) withlimT→∞1T∫ T0Ri(t, u) du = 0,although for this system we consider the periodic form of the Fredholm alternative12pi∫ 2pi0Ri(t, T ) dT = 0. (2.24)Both the general and periodic form of the Fredholm alternative have been well studiedand a more theoretic approach to the periodic version is discussed in Bensoussan’s Asymp-totic analysis for periodic structures [3]. From (2.21), we learn the leading order term isonly dependent on the slow time, x0 = x0(t). Applying the Fredholm alternative (2.24) to(2.22) gives an equation for the slow behavior which then implies an equation for the fastbehavior with0 =12pi∫ 2pi0(−x0t(t)− µ− 2x0(t) + x0(t)2 +A sin(T )) dT,x0t =− µ− 2x0 + x20, x1T = A sin(T ).(2.25)Solving for the equilibrium solution of (2.25) leads to the leading order solution, x0, andalso allows us to partially solve for the first correction term x1 withx0 =1−√1 + µ,x1(t, T ) =v1(t)−A cos(T ).Repeating this procedure in (2.23), as shown in Appendix B, results in the expansion (2.20)written in the original variablesx ∼ 1−√1 + µ− Ω−1A cos(Ωt) +O(Ω−2). (2.26)Once again, the explicit outer solution (2.26) performs well for x away from the axis x = 0,we search for when the assumptions of the asymptotic series fail indicating where an inneranalysis is needed. This is when x0 ∼ x1 which occurs for µ ∼ O(Ω−1).We consider a general scaling in the form of x = Ω−αy and µ = Ω−βm where α > 0and β > 0 allow for an inner equation to be found. Applying these local variables to (2.1)results iny˙ = −Ωα−βm+ 2|y|−Ω−αy|y|+ΩαA sin(Ωt). (2.27)In the local system (2.27), we find similar behavior occurring across multiple scales andthus we are able to use the same time scales from the outer analysis, t = t and T = Ωt.Then we have y(t) = y(t, T ) and hence a similar multiple scales argument in (2.27) leads toyT + Ω−1yt = −Ωα−β−1m+ Ω−12|y|−Ω−α−1y|y|+Ωα−1A sin(T ). (2.28)172.3. High Frequency Oscillatory ForcingWith a standard balancing argument between the leading order terms in (2.28) yT andΩα−1A sin(T ), we see that α = 1. We also want to see the terms Ωα−β−1m balance withΩ−12|y|, which gives us that β = 1 as well. This results in the inner equationyT + Ω−1yt = Ω−1 (−m+ 2|y|)− Ω−2y|y|+A sin(T ). (2.29)Similarly to the outer equation, we approximate the solution with an asymptotic expansionin terms of Ω−1y(t, T ) ∼ y0(t, T ) + Ω−1y1(t, T ) +O(Ω−2). (2.30)Substituting the expansion (2.30) into the inner equation (2.29) we findy0T + Ω−1y0t + Ω−1y1T + . . . = Ω−1(−m+ 2|y0 + Ω−1y1 + . . . |) +A sin(T )+ Ω−2(y0 + Ω−1y1 + . . .)|y0 + Ω−1y1 + . . . |.Here we then find the following system of equations at each order of ΩO(1) : y0T = A sin(T ), (2.31)O(Ω−1) : y1T + y0t = −m+ 2|y0|. (2.32)Solving the leading order equation (2.31) gives that the leading order term has the form,y0(t, T ) = v0(t)−A cos(T ). Applying the Fredholm alternative (2.24) to (2.32) leads tov0t(t) = −m+1pi∫ 2pi0|v0(t)−A cos(T )| dT. (2.33)In this setting, we must consider two cases for v0(t) that determine the nature of thisintegrand. Case I: if v0(t) is large enough to keep y0 from ever changing sign and Case II:if v0(t) is too small and y0 crosses the x = 0 axis. In figure 2.3 we show the range of eachcase, the region on the right is following under case I, the green dotted vertical line definingthe parameter range between the cases, the middle region for case II and the blue verticalline giving the bifurcation, µosc, which is determined below.182.3. High Frequency Oscillatory ForcingFigure 2.3: The parameter ranges for each case are shown here with A = 2. For reference,the original bifurcation diagram is overlayed.2.3.1 Case I: v0(t) ≤ −|A|We call this the ’below axis’ case as the solution stays far from the axis x = 0 for most of theoscillation and thus the behavior is not influenced by the non-smooth dynamics. We do notexpect to see the bifurcation occur under these conditions but instead we find the parameterrange for each of these cases. Here the integral in equation (2.33) is straightforward toevaluate as v0(t) is a constant with respect to the fast time T , thus we find the innerequation and equilibriumv0t = −m− 2v0, v0 = −m2.This gives the leading order equilibrium solution with oscillations of the local variables forthis case which we write in terms of the original variablesy(t, T ) ∼− m2−A cos(T ) +O(Ω−1),x(t) ∼− µ2− Ω−1A cos(Ωt) +O(Ω−2).(2.34)The condition v0(t) ≤ −|A| combined with the equilibrium allows us to establish when(2.34) holdsµ ≥ 2|A|Ω. (2.35)Following the equilibrium to (2.35) leads us to case II where the oscillations cross the axisand the assumptions of this case no longer hold.192.3. High Frequency Oscillatory Forcing2.3.2 Case II: |v0(t)|< |A|We call this the ’crossing’ case; here the equilibrium is small enough that the oscillationscan now push the solution above the axis. Under these conditions, the solution spends timenear the axis x = 0 and thus experiences non-smooth influence. As the crossing continues,the non-smooth behavior drives the solution to gradually grow. Therefore we expect to findthe bifurcation here. From (2.35), we have a range of µ when this case applies, µ < 2|A|Ω . Itis important to note that the integrand in (2.33) is non-trivial when |v0(t)|< |A|. In orderto deal with the sign changing inside the integral, we break the integration into regionsbased on the sign. Recall that we are searching for equilibrium behavior, and so we maymake the assumption that we are dealing with a fixed value of v0 such that |v0|≤ |A|. Infigure 2.4 we observe the function that we are integrating.Figure 2.4: The non-smooth function |y0(T )|= |v0−A cos(T )| that we integrate is shown asa solid red line. We also show an example of v0 as a horizontal blue dotted line. Here thevalue of |v0|≤ |A|, which causes kinks to appear at the roots of |y0|: T1 and T2 respectively.These are the vertical black dashed dotted lines.From figure 2.4, the roots of the integrand areT1 = arccos(v0/A), T2 = 2pi − arccos(v0/A).Here we notice 0 < T1 < T2 < 2pi and that the sign of the integrand stays the same on eachinterval. We only assume that the first interval [0, T1] observes a solution while the centeris still negative, thus the integrand will also be negative. From this, the integral in (2.33)202.3. High Frequency Oscillatory Forcingis computed as∫ 2pi0|v0 −A cos(T )| dT =−∫ T10(v0 −A cos(T )) dT+∫ T2T1(v0 −A cos(T )) dT −∫ 2piT2(v0 −A cos(T )) dT.(2.36)Evaluating (2.36) and using a trigonometric identity, sin(arccos(x)) =√1− x2, we find theintegral to be ∫ 2pi0|v0 −A cos(T )| dT = 2pi(arcsin(v0/A)v0 +√A2 − v20).Notice that our argument above is simple for fixed v0, but we have used a multiple scalesapproach for our fast time with t T implying that v0(t) is approximately fixed overT ∈ [0, 2pi]. This holds true due to having a high frequency Ω and otherwise would not bea valid approximation. Thus we can evaluate (2.33) to find the inner equationv0t = −m+4pi(arcsin(v0/A)v0 +√A2 − v20). (2.37)In its current form, (2.37) prevents v0(t) to be found analytically, so we use a quadraticTaylor approximation to be able to solve this equation explicitly. This then givesv0t ≈ −m+4|A|pi+2pi|A|v20, (2.38)which has the following equilibrium with positive constant Cv0 = −C√m− 4|A|pi. (2.39)Thus we have the leading order inner equilibrium (2.39) and writing this in the originalvariables givesy ∼− C√m− 4|A|pi−A cos(T ) +O(Ω−1),x(t) ∼− C√Ω(µ− 4|A|piΩ)− Ω−1A cos(Ωt) +O(Ω−2).(2.40)It then is clear that the bifurcation, µosc, occurs when (2.40) fails to be real valued. Thuswe find µosc to take the formµosc =4|A|piΩ. (2.41)From the result (2.41), we gather that the oscillatory forcing in the system causesthe bifurcation to occur sooner, µosc > µns, and this is controlled by the size of A andΩ. Heuristically, the model experiences the non-smooth behavior sooner in µ with theoscillations, but we can see as Ω → ∞ then µosc → µns. This effect is contrary to the212.3. High Frequency Oscillatory Forcingslow variation where the solution experienced a delayed tipping, µslow < µns. Although ourassumption that Ω 1 must be met for our analysis to hold, we may still see a medium-range Ω. From figure 2.5 (d) we may have Ω = 5 and get reasonable approximationswith µosc. This means that our approximation µosc improves with increasing Ω, but itstill provides a reasonable approximation for bifurcations away from µns for a range of Ω.Advanced bifurcation also indicates that the region of bi-stability is shrunk with oscillatoryforcing and thus can be used to eliminate the region entirely with A and Ω chosen properly,effectively destroying any hysteresis. We compare our estimate to numerical results forvarying sizes of Ω−1.(a) (b)(c) (d)Figure 2.5: In (a) the numerical time series solutions to (2.1) are given from bottom to topwith µ = {.8, .33, .15} in case I, case II and µ < µosc in (2.41) respectively with A = 2,Ω = 10 and = 0. In (b) we show the time series on the bifurcation diagram. In (c), azoom in closer to the non-smooth bifurcation of (b), where the dotted vertical lines dictatethe region between case I and case II (green) as well as the bifurcation estimate (blue)respectively. In (d) a range of Ω−1 and the corresponding numerical bifurcations (redstars) are compared to our estimate of the bifurcations (black solid line). We consider thebifurcation criterion to be when the numerical solution has passed x > .5.222.3. High Frequency Oscillatory ForcingIn figure 2.5 an example is given of the effect oscillatory forcing has on the solution givena choice of A and Ω with (a), (b) and (c), but (d) shows the bifurcation approximationacross a range of Ω−1. There is an allowed range of Ω from our assumption of Ω 1and in this region we see agreement between the numerical results and our approximations.The concavity is well represented and the behavior as Ω−1 → 0 converges to the staticbifurcation. Thus we expect that our methodology is valuable for the Stommel model.2.3.3 StabilityOnce more, the outer solution (2.26) is stable from the static model in section 2.1. Inthis section, we have two regions of interest and establish their stabilities agree with ouranalysis. Each region has a particular version of the same inner equation dictating thesolution’s behavior, namelyv0t = −m+1pi∫ 2pi0|v0 −A cos(T )| dT. (2.42)Case I: v0(t) < −|A|In this region we did not find any bifurcation behavior and (2.42) simplifies to the innerequation with equilibrium z0 as followsv0t = −m− 2v0 = f(v0), z0 = −m2.Similarly to section 2.2, we have a fixed parameter equation that we have shown to causeperturbations to decay exponentially and hence we find the equilibrium to be hyperbolicand asymptotically stable in the paramter range found in the analysis (2.35)µ ≥ 2|A|Ω.Case II: |v0(t)|< |A|For this region we found the bifurcation and hence we should expect to lose stability here.The Taylor approximation (2.38) for the inner equation with equilibrium z0 isv0t = −m+4|A|pi+2pi|A|v20, z0 = −C√m− 4|A|pi. (2.43)We consider a simple linear perturbation of (2.43), v0(t) = z0 + u(t) with ‖u(t)‖ 1.Applying the standard Taylor expansion to determine the equation for the perturbations,we findv0t =f(z0) + fv0(z0)(v0 − z0) +O(‖v0 − z0‖2),ut =− 2√m− 4|A|piu.(2.44)The sign of (2.44) gives that the perturbations decay exponentially and hence the equi-librium is hyperbolic and asymptotically stable as long as m > 4|A|pi or equivalently µ >4|A|piΩ .232.4. Slowly Varying and Oscillatory ForcingWe find that once µ reaches this value in (2.41), then the stability of (2.44) is non-hyperbolic.When the stability switches like this, we expect a bifurcation. Thus we have further evidenceto support that (2.41) is the oscillatory bifurcation we seek.2.4 Slowly Varying and Oscillatory ForcingNow that we have established an approach for each feature of the model individually, wecombine them in the full one component model (2.1) where 1 and A ∼ O(1). Dueto the slow variation in µ, we do not see a bifurcation occur under these conditions butrather a tipping point. Hence we must find the behavior of the solution and search forwhen a rapid transition towards the upper branch occurs. Since the high frequency couldbe approximated by a power of slow variation, we choose to relate these mechanisms witha generic polynomial Ω = −λ for a parameter λ > 0. With a general λ, we classify regionsof behavior by ranges of λ and are able to determine where mixed behavior occurs orwhen one of the mechanisms becomes dominant. With both mechanisms in effect, we againchoose to use a multiple scales approach to capture both slow behavior and fast oscillations.Although now, we truly have slow behavior, the slowly varying parameter µ(t), as well asfast behavior, the rapid oscillations sin(Ωt). The choice in time scales is then τ = t andT = −λt, which leads to the systemxT + λ+1xτ = λ(−µ(τ) + 2|x|−x|x|+A sin(T )),µτ =− 1.Once again, we assume initial conditions satisfying x < 0 and starting far enough awayfrom x = 0, before any crossing occurs, to find the outer solution. Thus we have the systemxT + λ+1xτ = λ(−µ(τ)− 2x+ x2 +A sin(T )),µτ (τ) =− 1.(2.45)We perform an asymptotic expansion in terms of the small quantity λ, where we note thatthis is the same as Ω−1,x(τ, T ) ∼ x0(τ, T ) + λx1(τ, T ) +O(1+λ, 2λ). (2.46)Introducing the expansion (2.46) into the outer multi-scaled equation (2.45) givesx0T +λ+1x0τ+λx1T +. . . = λ(−µ(τ)−2x0+x20+A sin(T ))+2λ(−2x1+x1x0)+. . . (2.47)Here we separate (2.47) at each order of λ to find the following system of equationsO(1) : x0T = 0, (2.48)O(λ) : x1T = −µ(τ)− 2x0 + x20 +A sin(T ), (2.49)O(2λ) : x2T + 1−λx0τ = −2x1 + 2x0x1. (2.50)242.4. Slowly Varying and Oscillatory ForcingDepending on the value of λ, O(λ+1) may be the next order before O(2λ), althoughconsidering either produce the same equation at their respective order and hence our choicein λ does not change the calculations up to this correction term. Thus for the outer solution,we consider the system with O(2λ). Each equation gives the behavior of each order ofthe solution; (2.48) indicates that the leading order term is only slow time dependent,x0 = x0(τ). In Appendix B we apply the Fredholm alternative (2.24) to (2.49) and (2.50)to find the first few terms of the expansion in (2.46) explicitly. The resulting solution isx(t) ∼ 1−√1 + µ(t)− 4(1 + µ(t))− λA cos(Ωt) +O(1+λ, 2λ). (2.51)In the outer solution (2.51), we consider when the terms violate the assumptions of theexpansion to find where we need to use an inner equation. This happens either whenx0 ∼ O() or when x0 ∼ O(λ) which is when µ ∼ O() or µ ∼ O(λ) respectively anddepends on the value of λ.To find an inner equation we use a general scaling for both x and µ given the ambiguityof the choice in µ withx(t) = αy(t), µ(t) = βm(t), (2.52)where α > 0 and β > 0 allow for inner equations to be found. Applying the local variables(2.52) to the canonical equation (2.1) givesαy˙ =− βm(t) + α2|y|−2αy|y|+A sin(−λt),m˙ =− 1−β.(2.53)From (2.53) we find the fast time still appears but the slow time has multiple choicesdepending on λ. For convenience we choose to take a multiple scales approach with scalest and T = −λt in (2.53) to findα−λyT + αyt =− βm(t) + α2|y|−2αy|y|+A sin(T ),mt =− 1−β.(2.54)To determine the correct scalings in (2.52), we balance the leading order terms on bothsides of (2.54) α−λyT and A sin(T ), which gives us that α = λ. This suggests that theoscillatory term persists in the inner asymptotic expansion of (2.1) regardless of the choicein λ.We now consider the same scales t and T = −λt on the canonical system (2.1)xT + λxt =− λ+βm(t) + λ2|x|−λx|x|+λA sin(T ),mt =− 1−β.(2.55)Here we use the expansionx(t, T ) = λy0(t, T ) + . . . ,where the next terms of this expansion depend on whether λ ≤ 1 or λ > 1. We considerthese ranges in case I and case II respectively.252.4. Slowly Varying and Oscillatory Forcing2.4.1 Case I: λ ≤ 1We call this the ’mixed effects’ case a both slow variation and oscillatory forcing causingnoticeable effects on the solution for this range of λ. Hence we consider the expansionx(t, T ) ∼ λy0(t, T ) + qy1(t, T ) + . . . , (2.56)with the next term q > λ to be consistent with the scale analysis above that determinedinner behavior to start at O(λ). Substituting (2.56) into (2.55) givesy0T + λy0t + q−λy1T + qy1t + . . . = − βm(t) + λ2|y0 + q−λy1 + . . . |+A sin(T )+ 2λ(y0 + q−λy1 + . . .)|y0 + q−λy1 + . . . |.Separating by distinct orders of then gives the following equations at each orderO(1) : y0T = A sin(T ), (2.57)O(λ) : q−2λy1T + y0t = −β−λm(s) + 2|y0|. (2.58)In (2.58) we find that the appropriate next term in the expansion (2.56) is with q = 2λ.This choice in q keeps the equations balanced but q implies that λ > 12 for an expansionto be found. Otherwise, the quadratic terms must be included and we no longer find localequations. This indicates that the range of λ ≤ 12 behaves differently. We discuss this furtherin chapter 4. There is also the choice between β = λ or β = 1 and each has a particularappeal. With β = λ, the form of (2.58) is simple, but the equation for the slow variation ismt = −1−λ. This then suggests a slower time scale to approach the problem. We insteadchoose to allow β = 1 for convenience and track a small coefficient on m(t) in exchangefor keeping the same time scale with mt = −1. This is valid as long as we are trackingsmall coefficients and not large ones as this would suggest we used an incorrect scaling, butboth of these choices lead to the same conclusion. Using (2.57) gives the appropriate form,y0(t, T ) = v0(t)−A cos(T ). We then apply the Fredholm alternative (2.24) to (2.58) whichgives a similar equation to the integral (2.33) in section 2.3 withv0t = −1−λm(t) +1pi∫ 2pi0|v0(t)−A cos(T )| dT. (2.59)The approach developed in section 2.3 is applied here to (2.59), where we separate thebehavior of the integral based on the relative size of v0(t) to A. We have the followingsituations, sub-case I: v0(t) ≤ −|A| and sub-case II: |v0(t)|< |A|.Sub-Case I: v0(t) ≤ −|A|Once more, we call this the ’below axis’ sub-case and we do not expect tipping to occurunder these conditions since the solution is entirely negative and far from the axis x = 0for most of the oscillation. Under these conditions, (2.59) gives the simple inner equationv0t = −1−λm(t)− 2v0. (2.60)262.4. Slowly Varying and Oscillatory ForcingSolving (2.60) can be done under our assumptions much like in subsection 2.3.1 but insteadwe focus on the pseudo-equilibrium. This choice results in finding the effective parameterrange for µ which distinguishes these sub-cases and helps to determine when the solu-tion enters sub-case II. Since m(t) is allowed to vary, this must be thought of more as apseudo-equilibrium and we are only interested in when the pseudo-equilibrium violates theassumptions of this case. Finding the pseudo-equilibrium of (2.60) givesv0(t) = −1−λm(t)2.Using the condition v0(t) ≤ −|A| gives that m(t) ≥ λ−12|A|. Writing this result in originalvariables gives us the parameter rangeµ(t) ≥ 2|A|Ω, (2.61)for sub-case I which agrees with the range from (2.35) in section 2.3. Following the pseudo-equilibrium to the boundary (2.61), we eventually reach sub-case II where we see the oscil-lations crossing the axis.Sub-Case II: |v0(t)|< |A|Again, we call this the ’crossing’ sub-case. Here the behavior of the solution dependsstrongly on the sign of the solution similarly to section 2.3. We seek the relationship betweenslow variation and oscillatory forcing on the tipping point. As the pseudo-equilibrium getscloser to the x = 0 axis, the solution spends more time above this axis and more complicatedcontributions from the sign changing appear. We expect tipping to happen under theseconditions.The methodology of solving the integral in (2.59) holds identically to that of subsec-tion 2.3.2. Here, we have a slow time function v0(t) that is approximately fixed withrespect to the fast time T under the multiple scales approach. Thus we evaluate the in-tegral by separating the sign of the integrand with the values T1 = arccos(v0/A) andT2 = 2pi − arccos(v0/A) to findv0t = −1−λm(t) +4pi(arcsin(v0/A)v0 +√A2 − v20). (2.62)We then choose to find an explicit analytic expression by approximating (2.62) with aquadratic Taylor expansion. This givesv0t =− 1−λm(t) +4|A|pi+2pi|A|v20,mt =− 1.(2.63)With (2.63) in terms of slow time, it restricts any analytical approaches that link the effectsof the varying parameter. Instead we take the same approach from [7] and switch thedifferentiation onto the slow varying parameter m withv0m = 1−λm− 4|A|pi− 2pi|A|v20. (2.64)272.4. Slowly Varying and Oscillatory ForcingIt is here where we take advantage of the form of (2.64) with the form from (1.3) to solve,resulting inv0(m) ∼ (1−λ)/3(pi|A|2)2/3 Ai′(2(λ−1)/3 ( 2pi|A|)1/3 (1−λm− 4|A|pi ))Ai(2(λ−1)/3(2pi|A|)1/3(1−λm− 4|A|pi )) .With the solution to (2.56) we rewrite back into the original variablesy0(t, T ) ∼CAi′(2(λ−1)/3(2pi|A|)1/3(1−λm(t)− 4|A|pi ))Ai(2(λ−1)/3(2pi|A|)1/3(m(t)− 4|A|pi )) − λA cos(T ) + . . . ,x(t) ∼CAi′((Ω2)1/3 ( 2pi|A|)1/3(µ(t)− 4|A|piΩ ))Ai((Ω2)1/3 ( 2pi|A|)1/3(µ(t)− 4|A|piΩ )) − λA cos(Ωt) + . . .(2.65)Given the inner solution (2.65), we search for the singularity of this solution in order toidentify tipping. Recall from (1.4) that the singularity relates to the first root of the Airyequation. Here we find the singularity µmixed to beµmixed =(2Ω)1/3(pi|A|2)1/3(−2.33811 . . .) + 4|A|piΩ. (2.66)The value µmixed which causes this singularity is our tipping point. We rewrite (2.66) toemphasize the contributions from the slow variation of the parameter and the oscillatoryforcingµmixed =(pi|A|2Ω)1/3µsmooth + µosc, (2.67)with µsmooth = 2/3 (−2.33811 . . .), similarly to the smooth problem from [24], and µosc from(2.41) respectively.The resulting tipping approximation (2.67) indicates that the size of the amplitude Adetermines whether the tipping occurs early or late relative to the bifurcation. Naturallywe see a larger amplitude cause more contribution from the oscillations and hence an earliertipping. On the other hand, larger values in cause this tipping to occur later. So theseeffects have opposite pulls on the tipping and can effectively cancel one another out underproper conditions. It would even be possible to break the hysteresis cycle by eliminatingthe region of bi-stability in this model with sufficiently large amplitude and small . Thetipping point holds for any λ ∈ (12 , 1] and we see different behavior for larger λ.2.4.2 Case II: λ > 1We call this the ’slowly varying dominant’ case as this is when we see that the oscillationscontribute less than the slow variation. For this range of λ the scaling for µ is simple,282.4. Slowly Varying and Oscillatory Forcingµ = m. Thus we expect to see integer powers in the leading order along with powers of λso we choose the expansionx(t, T ) ∼ λy0(t, T ) + y1(t, T ) + qy2(t, T ) + . . . , (2.68)where q > λ to allow for consistency with the scale analysis but not necessarily the samevalue as in case I. Substituting (2.68) into (2.55) givesy0T + λ+1y0t + λy1T + qy2T + . . . =− λ+1m(t) + λ+12|y0 + λ−1y1 + . . . |+ λ+2(y0 + λ−1y1 + . . .)|y0 + λ−1y1 + . . . |+ λA sin(T )Here we separate out each order of to find the equations at each orderO() : y0T = 0, (2.69)O(λ) : y1T = A sin(T ), (2.70)O(λ+1) : q−λ−1y2T + y0t = −m(t) + 2|y0 + λ−1y1|. (2.71)We learn in (2.71) that q = λ + 1 keeps the terms balanced. From (2.69) we find thatthe dominant behavior for this case is only slow time dependent, y0 = y0(t) and from(2.70) that the oscillatory behavior occurs in y1 with y1(t, T ) = v1(t) − A cos(T ). Sincewe have y1 as a correction to y0, we may absorb the slow behavior into y0. Thus we treaty0(t) = y0(t) + λ+1v1(t) ≈ y0(t). Applying Fredholm to (2.71) givesy0t =−m(t) +1pi∫ 2pi0|y0(t)− λ−1A cos(T )| dT. (2.72)With λ ≈ 1, we see nearly identical behavior in (2.72) as that of what we explored insubsection 2.4.1. As long as the amplitude of oscillations inside the integral are λ−1A ∼O(1), then this integral is similar to the integral in Case I (2.59). To see this, we followthe same approach as to integrate (2.72) with T1 = arccos(y0/λ−1A)and T2 = 2pi −arccos(y0/λ−1A)which givesy0t = −m(t) +4pi(arcsin(y0/λ−1A)y0 +√(λ−1A)2 − y20). (2.73)Here we apply the same quadratic Taylor approximation to (2.73) to findy0t =−m(t) + λ−12|A|pi+ 1−λ2pi|A|y20,m =− 1.(2.74)We again use the result from (1.4) to find the tipping, which we then write into originalvariables292.4. Slowly Varying and Oscillatory Forcingmmixed =(λ−1)/3(pi|A|2)1/3(−2.33811 . . .) + λ−1 4|A|pi,µmixed =(pi|A|2Ω)1/3µsmooth + µosc.Thus we conclude that there is a natural transition into case II from case I with almostthe same behavior and identical tipping as in (2.67). As λ continues to grow, the amplitudeof oscillation in (2.72) decays and the contribution from the oscillations weaken. This allowsus to say that the integral is approachingy0t = −m(t) + 2|y0|. (2.75)With (2.75) taking the same form as in section 2.2, this allows us to use the results thereto find the solution. We write this in terms of the original variablesy0(t, T ) ∼Ce−2m(t) + m(t)2− 1/4 + λA cos(T ),x(t) ∼Ce−2µ(t)/ + µ(t)2− λA cos(Ωt) +O(2λ).(2.76)This then leads to the same tipping as in the slowly varying model withµslow =12 log .Thus we find that in this case and with λ near 1, the same tipping point µmixed in (2.67) fromsubsection 2.4.1 is found. For large λ, the oscillations have less of an impact and the solutiontips entirely like µslow in (2.14) from section 2.2. For convenience, this is summarized in thefollowing table.One Component Tipping Points > 0 and A = 0: µslow = ln()/2 = 0 and A 6= 0 with Ω 1: µosc = 4|A|piΩ > 0, A 6= 0 and λ ≤ 1: µmixed =(pi|A|2Ω)1/3µsmooth + µosc > 0, A 6= 0 and λ > 1 with λ ≈ 1: µmixed =(pi|A|2Ω)1/3µsmooth + µosc > 0, A 6= 0 and λ > 1: µslow = ln()/2Table 2.1: Overview of the tipping points of the one component model for each mechanismand case.In figure 2.6, we see an example of the numerical solution to the canonical system (2.1)with slow variation and oscillatory forcing. This example has a tipping point occurring incase I due to λ ∈ (12 , 1] allowing the slow variation and oscillatory forcing to produce amixed effect on the tipping point. Although we see noticeable contributions from the slowvarying parameter the tipping point still occurs near the oscillatory bifurcation. This tells302.4. Slowly Varying and Oscillatory Forcingus that for these choices in the values the strongest effect is the oscillatory forcing. It ispossible to find values of , A and λ that cause the non-smooth tipping to occur at the sameplace as the smooth bifurcation. This would eliminate the region of bi-stability and destroythe hysteresis curve entirely for this model.(a) (b)Figure 2.6: On the left, one can see the bifurcation diagram for the canonical system (2.1)with the numerical solution (black dotted line). On the right, a zoom in around the non-smooth bifurcation. The dotted vertical line is the tipping point µmixed (2.66) (blue). Thevertical line (black) is the when the numerical solution has passed the tipping criterionx > .5. The parameter values are = .05, λ = .8 and A = 4.In figure 2.7, we see an example of λ falling into case II yet close enough to 1 that we seemixed behavior in the tipping. Here the slow variation is dominant and the oscillations areonly noticeable in the zoom in. The green dotted line is the tipping point approximation(2.14) from section 2.2 and µmixed is still approximating the tipping point well.312.4. Slowly Varying and Oscillatory Forcing(a) (b)Figure 2.7: On the left, one can see the bifurcation diagram for the canonical system (2.1)with the numerical solution (black dotted line). On the right, a zoom in about the non-smooth bifurcation. The dotted vertical lines are the tipping point µmixed (2.66) (blue)and slowly varying tipping µslow (2.14) (green). The vertical line (black) is the when thenumerical solution has passed the tipping criterion x > .5. The parameter values are = .05,λ = 1.05 and A = 4.In figure 2.8, we see an example of λ falling into case II but large enough that wesee almost entirely slow behavior in the tipping. Even upon closer inspection it is hardlynoticeable that oscillations are present in the model. The green dotted line is the tippingapproximation (2.14) from section 2.2, and it is clear that this is a better approximationthan the mixed tipping point.322.4. Slowly Varying and Oscillatory Forcing(a) (b)Figure 2.8: On the left, one can see the bifurcation diagram for the canonical system (2.1)with the numerical solution(black dotted line). On the right, a zoom in around the non-smooth bifurcation. The dotted vertical lines are the tipping point µmixed (2.66) (blue)and slowly varying tipping µslow (2.14) (green). The vertical line (black) is the when thenumerical solution has passed the tipping criterion x > .5. The parameter values are = .05,λ = 1.6 and A = 4.In figure 2.9 we compare the tipping between case I and case II with the numericaltipping. For smaller λ, the frequency Ω gets smaller and the case I tipping becomes morepredominant. For the analysis performed in this section, Ω 1 and for λ ≤ 12 we haveΩ ∼ O(1). We do not consider a low frequency corresponding to λ ≤ 12 in this thesis. Thelarger λ becomes, the less effect we see due to the oscillatory forcing until it is negligiblefor some λ > 1. This is also seen in the asymptotic solution for each case, (2.51), (2.65),and (2.76), where the oscillatory component of the term has a λ coefficient and shrinks theeffects as λ grows.332.4. Slowly Varying and Oscillatory ForcingFigure 2.9: An example of numerical tipping (red stars) as the numerical solution to (2.1)passes the tipping criterion x = .5 for the last time. The parameter values are = .01 andA = 4. The lines are the case I tipping estimate (2.66) (black solid line) and the case IItipping estimate (2.14) (blue dotted line).The performance of our estimates are seen in figure 2.10. For case I tipping, the rangeof appropriate is highly dependent on the choice in λ. Often, the range is quite small foraccurate estimates, but with this in mind both case approximations have good performanceover a reasonable range of .(a) λ = .8 (b) λ = 1.3Figure 2.10: The numerical tipping (red stars) follows the appropriate case depending on λ.The case I tipping estimate (black solid line) and the case II tipping estimate (blue dottedline) are shown.342.4. Slowly Varying and Oscillatory Forcing2.4.3 StabilitySimilarly to section 2.3, there are two ranges for λ that govern the stability of solutionsfound in our analysis, namely λ ≤ 1 and λ > 1.Case I: λ ≤ 1Recall that for this case, we must have λ ∈ (12 , 1] and for this range, we found the innerequationv0t = −1−λm(t) +1pi∫ 2pi0|v0(t)−A cos(T )| dT = f(t, v0). (2.77)As in our previous analysis, we consider two regions of the solution v0(t) for this integral,in subsection 2.4.1 with sub-case I: v0(t) ≤ −|A| and in subsection 2.4.2 with sub-case II:|v0(t)|≤ |A| where each these sub-cases deal with the respective size of v0(t) to the amplitudeof the effective oscillations.Sub-Case I: v0(t) ≤ −|A|Recall from the analysis that equation (2.77) simplifies in this region of v0(t) and has thefollowing inner equation and pseudo-equilibrium z0(t)v0t = −1−λm(t)− 2v0 = f(t, v0), z0(t) = −1−λm(t)2. (2.78)As we saw in section 2.2, special treatment of the pseudo-equilibrium stability analysis isneeded with linear perturbations v0(t) = z0(t)+u(t), where ‖u(t)‖ 1 and z0t = −1−λmt2 =1−λ2 . The resulting Taylor expansion is thusv0t =f(t, z0) + fv0(t, v0)(v0(t)− z0(t)) +O(‖v0(t)− z0(t)‖2),ut =− 1−λ2− 2u.This leads to the conclusion that equation (2.78) causes perturbations to decay exponen-tially to just below the pseudo-equilibrium. Hence we find the pseudo-equilibrium to be ahyperbolic and asymptotically attracting solution.Sub-Case II: v0(t) ≤ |A|With the Taylor approximation from the analysis (2.63), we have the following inner equa-tion and pseudo-equilibrium z0(t)v0t = −1−λm(t) +4|A|pi+2pi|A|v20 = f(t, v0), z0(t) = −C√1−λm(t)− 4|A|pi. (2.79)352.4. Slowly Varying and Oscillatory ForcingWe consider simple linear perturbations to this pseudo-equilibrium (2.79), v0(t) = z0(t) +u(t) with ‖u(t)‖ 1. Treating the pseudo-equilibrium carefully, we find that the slowlyvarying component of the equilibrium contributes to the derivative. Thus we havev0t =z0t (t) + ut,z0t (t) =1−λ2C√1−λm(t)− 4|A|pi1−λm(t) > 4|A|pi ,0 1−λm(t) = 4|A|pi .(2.80)Now applying a Taylor expansion, we find the following behavior of perturbationsv0t =f(t, z0) + fv0(t, z0)(v0 − z0(t)) +O(‖v0(t)− z0(t)‖),ut =−1−λ2C√1−λm(t)− 4|A|pi− 2√1−λm(t)− 4|A|pi u 1−λm(t) > 4|A|pi ,0 1−λm(t) = 4|A|pi .(2.81)From (2.81), we find that the perturbations decay to a fixed negative quantity. Thisindicates, much like in section 2.2, that there is an attracting solution below the pseudo-equilibrium. The negative sign describes exponential decay and hence this equilibrium ishyperbolic and asymptotically stable for 1−λm(t) > 4|A|pi or µ(t) >4|A|piΩ . For 1−λm(t) = 4|A|pior µ(t) = µosc, the stability of (2.81) suddenly becomes non-hyperbolic. This tells us thatwe lose stability at the oscillatory bifurcation but the tipping point occurs afterwards, whichagrees with the conclusion in the tipping approximation from (2.67).Case II: λ > 1From the analysis, we discovered that as long as λ−1A ∼ O(1), then we have a similarbehavior in the tipping point. With the Taylor approximation from the analysis (2.72), theinner equation and pseudo-equilibrium z0(t) arey0 =−m(t) + λ−1 2|A|pi+ 1−λ2pi|A|y20 = f(t, y),z0(t) =− λ−1C√m(t)− λ−1 4|A|pi.(2.82)Similarly to Case I, we consider simple linear perturbations to this pseudo-equilibrium(2.85), y0(t) = z0(t) + u(t) with ‖u(t)‖ 1. Treating the pseudo-equilibrium carefully,we find that the slowly varying component of the equilibrium contributes to the derivative.Thus we havey0t =z0t (t) + ut,z0t (t) =λ−12C√m(t)−λ−1 4|A|pim(t) > λ−1 4|A|pi ,0 m(t) = λ−1 4|A|pi .(2.83)362.4. Slowly Varying and Oscillatory ForcingNow applying a Taylor expansion, we find the following behavior of perturbationsy0t =f(t, z0) + fy0(t, z0)(y0 − z0(t)) +O(‖y0 − z0(t)‖2),ut =−λ−12C√m(t)−λ−1 4|A|pi− 2√m(t)− λ−1 4|A|pi u m(t) > λ−1 4|A|pi ,0 m(t) = λ−1 4|A|pi .(2.84)The conclusions from case I still apply to (2.84) and thus we still have an attracting solutionuntil µ(t) = µosc and expect to see tipping occurring after the oscillatory bifurcation whichis consistent with our tipping approximation for this case.On the other hand, for large λ the integral (2.72) approachesy0t = −m(t) + 2|y0|. (2.85)This is the same type of behavior from section 2.2, where we found that for m(t) ≥ 0 ourpseudo-equilibrium was attracting and for m(t) < 0 searching for the pseudo-equilibriumcaused a contradiction. Thus, we conclude that the tipping point occurs in the region ofm(t) < 0 which agrees with (2.14).37Chapter 3Two Component ModelWith the methods and approaches developed in chapter 2 for the one component model, wehave an expectation of the behavior of the two component Stommel model around the non-smooth bifurcation under similar conditions. Recall that we study the non-dimensionalizedmodel (1.6) and η2 is the control parameter linked to the flow rate. With the bifurcationstructure we explored in chapter 1, we consider a generalization of the Stommel modelV˙ = η1 − η2 + η3(T − V )− T − V |V |+A sin(Ωt),T˙ = η1 − T (1 + |V |) +B sin(Ωt),η˙2 = −V (0) = V 0, T (0) = T 0, η2(0) = η20,(3.1)with slow variation 1, high frequency Ω 1, amplitudes of oscillation A and B,and model parameters η1 and η3 as positive constants. We assume the initial conditionsV 0, T 0 and η02 > η2ns are along the lower equilibrium branch which puts emphasis onthe non-smooth behavior of the Stommel model with two additional features. First, weallow for slow variation in the bifurcation parameter which has been shown to be realisticsince η2 is related to the freshwater flux and therefore not a fixed parameter; the sameassumption is made in [19]. Second, we consider periodic forcing in the additive parameterη1 to account for oscillations in seasonal effects, annual effects as well as tidal currents.This same approach is taken in [16] to look into idealized forcing in the Atlantic MeridonalOverturning Circulation (AMOC). It should be noted that although oscillatory forcing ispresent, there is strong evidence to believe that stochastic forcing is also present in theAMOC [4] and [16] which is suggestive of the same being true for the THC. Thus weconsider periodic forcing in the Stommel model, but we discuss stochastic forcing in thefuture work section. To fully understand the effects of each component on the model, weconsider them individually before putting them together.For the remainder of this thesis, we make two assumptions: first that η3 < 1 whichgives the smooth bifurcation Vsmooth of (3.1) in the region V > 0. The value of η3 describesthe relative strength of the temperature relaxation to that of salinity, and it is frequentlyassumed that salinity has a slower relaxation, giving η3 < 1. Thus we restrict our attentionto η3 < 1 and the case of η3 > 1 follows similarly. The second assumption we make is thateven though we have a model in terms of V and T , the variable V is driving the dynamicsof the system as confirmed in the analysis below. This assumption means that we want tounderstand the non-smooth behavior in V where T follows in response to the effects of V .The evidence below shows that change in temperature respond to change in salinity. Thisassumption reduces the model, expressing the behavior of T in terms of V to find equationsin only one variable.383.1. Slowly Varying Bifurcation Parameter3.1 Slowly Varying Bifurcation ParameterWe consider the slow variation mechanism to understand the effects on the Stommel model(3.1) with 1 and A = B = 0, where the bifurcation parameter η2 slowly varies withoutoscillatory forcing. We expect to find a tipping point in the neighborhood of the aforemen-tioned non-smooth bifurcation at η2ns = η1η3. With the choice of η3 < 1, the lower branchwith V < 0 is the branch we focus on in order to approach the non-smooth behavior, thus(3.1) isV˙ = η1 − η2(t) + η3(T − V )− T + V 2,T˙ = η1 − T (1− V ),η˙2 = −.(3.2)From section 2.2, we learned that the one component model had a solution that displayedone type of behavior away from the non-smooth axis x = 0 and another type close to x = 0thus a local analysis gave insight into the tipping. Here we search for an outer solutionto (3.2) that helps us understand the behavior of the system away from the non-smoothV = 0. Since we have slow variation in η2, we choose to scale the system (3.2) with theslow time τ = tVτ =η1 − η2(τ) + η3(T − V )− T + V 2,Tτ = η1 − T (1− V ),η2τ = −1.(3.3)We use asymptotic expansions in terms of the small quantity to look for slowly varyingsolutions. Here we chooseV (τ) ∼V0(τ) + V1(τ) + 2V2(τ) + . . . ,T (τ) ∼T0(τ) + T1(τ) + 2T2(τ) + . . . ,(3.4)and substitute (3.4) into (3.3) to findV0τ + 2V1τ + . . . =η1−η2(τ) + η3(T0 − V0)− T0 + V 20+ (η3(T1 − V1)− T1 − 2V1V0) + . . .T0τ + 2T1τ + . . . =η1 − T0(1− V0) + (−T1(1− V0) + V1T0) + . . . .Separating at each order of then gives the following system of equationsO(1) :{0 = η1 − η2(τ) + η3(T0 − V0)− T0 + V 20 ,0 = η1 − T0(1− V0),(3.5)O() :{V0τ = η3(T1 − V1)− T1 + 2V1V0,T0τ = −T1(1− V0) + V1T0.(3.6)393.1. Slowly Varying Bifurcation ParameterWe then solve (3.5) simultaneously for the pseudo-equilibria noting that T0 is explicitly interms of V0T0(V0) =η11− V0 ,0 = η1 − η2(τ)− T0(V0) + η3(T0(V0)− V0) + V 20 .(3.7)With T0 and V0 found, we solve (3.6) for T1 explicitly in terms of V1T1(V1) =T0τ − T0(V0)V11− V0 ,V1 =−(1− V0)V0τ + (1− η3)T0τ (V0τ )(1− η3)T0(V0) + (η3 − 2V0)(1− V0) .(3.8)This gives the first few terms of the asymptotic expansion (3.4) with (3.7) and (3.8). Giventhese expressions, in the two component problem it is not immediately obvious when theouter solution breaks down. However, noting that for V0 → 0 the asymptotic expansion forV is no longer valid, we rescale the system to find an inner equation in this region.Analogously to section 2.2, this corresponds to the non-smooth bifurcation at (η2ns, Vns, Tns) =(η1η3, 0, η1), so we rescale (3.1) around these values. This results in the local variablesη2 =η2ns + n,V =X,T =η1 + Y.(3.9)Then we introduce these local variables (3.9) into the Stommel model (3.1) to find thefollowing inner systemX˙ = −n(t)− η3X − (1− η3)Y − X|X|,Y˙ = −η1|X|−Y − |X|Y,n˙ = −1.(3.10)Here we outline the influence of the parameters η1 and η3 on the behavior of a solution. Wealready determined in the introduction that η3 determines the orientation of the problem.By viewing (3.10) as a 2× 2 system of spatial variables, we find an interaction between theparameters η1 and η3 in the eigenvalues of this system. Linearizing then gives(X˙Y˙)=( −η3 −(1− η3)−η1sgn(X) −1)(XY)−(n(t)0). (3.11)A linearized stability analysis about the pseudo-equilibria similar to [21] is needed todetermine the stability of the lower pseudo-equilibria. This is performed in [5] and hence wenote that both stability of the lower branch and non-hyperbolic behavior at (η2ns, Vns, Tns)is observed. Thus we expect our tipping point to occur just after the static non-smooth403.1. Slowly Varying Bifurcation Parameterbifurcation η2ns. With this in mind, we consider V > 0 and with (3.11) we find the innersystemX˙ = −n(t)− η3X − (1− η3)Y − X2,Y˙ = −η1X − Y − XY,n˙ = −1.(3.12)Following the approach from section 2.2, we change the differentiation to be in terms of nto find (XnYn)=(η3 1− η3η1 1)(XY)+(n+ X2XY).Seeking a leading order solution in this region, we drop the order terms to give(XnYn)=(η3 1− η3η1 1)(XY)+(n0). (3.13)For the system (3.13), we find the eigenvaluesλ1,2 =η3 + 12± 12√(1 + η3)2 + 4(η1(1− η3)− η3). (3.14)The eigenvalues in (3.14) must be real as η3 < 1 guarantees the discriminant is alwayspositive. The type of stability also follows as one of the eigenvalues is positive and theother is negative. This causes the solution to be unstable, which confirms tipping to occurin the region V > 0. With real eigenvalues, the solution in the V > 0 region takes thefollowing exponential form with constants Ki,j being component j of the corresponding i-theigenvectorX(n) ∼K1,1eλ1n +K2,1eλ2n + C1n+ C2,Y (n) ∼K1,2eλ1n +K2,2eλ2n + C3n+ C4.(3.15)Translating both solutions in (3.10) back to our original variables we findV (t) ∼K1,1eλ1(η2(t)−η1η3)/ + K2,1eλ1(η2(t)−η1η3)/ +O(),T (t) ∼η1 + K1,2eλ1(η2(t)−η1η3)/ + K2,2eλ2(η2(t)−η1η3)/ +O().(3.16)With (3.16) admitting the solution in the region V > 0, we determine the system totip once one of these exponentials becomes large (i.e O(1/)). This causes the system toabruptly transition away from our lower branch towards the upper stable branch. Thetipping point η2slow is thenη2slow = mini{η2ns − ln()/λi}, i = 1, 2. (3.17)Thus we have the tipping for this problem with (3.17) and this has a similar form to thetipping point from section 2.2. As we found in (3.14), one of the eigenvalues is alwayspositive and thus η2slow < η2ns. This means the slow variation causes a delay in the rapidtransition from the lower branch to the upper branch and the solution spends more time413.1. Slowly Varying Bifurcation Parameteraround the lower branch. These effects shrink as → 0 until we return to the static problemwith = 0.(a) (b)Figure 3.1: In (a) the numerical solutions (black dotted and dash-dotted lines) to (3.1) aregiven with η1 = 4, η3 = .375, and = {.01, .04} respectively. In (b) a zoom is given closerto the non-smooth bifurcation region. The blue vertical lines are the predictions (3.17)against the black solid vertical lines which are the tipping points with the tipping criterionV > Vsmooth on the numerical solution.In figure 3.1 the numerical solution with a slowly varying η2 is given for two choices of in (a) and we zoom in around the non-smooth bifurcation in (b). Here we use the tippingcriterion to be whenever V > Vsmooth which is large enough that the numerical solution isstrongly moving towards the upper branch. The delay in moving towards the upper branchin the V solution causes a similar delay in T , as seen in figure 3.2. This is best seen withthe outer solution (3.4) as the correction terms are negative. Thus, when the tipping pointin V is reached, the solution for T has yet to reach the maximum value. Notice that afterthe tipping occurs, the numerical solution passes entirely over the unstable branch and evensome of the upper stable branch before following the upper branch closely.423.1. Slowly Varying Bifurcation Parameter(a) (b)Figure 3.2: In (a) we have the numerical solutions (black dotted and dash-dotted) over thestandard equilibrium plot for V vs. T . In (b) a zoom in closer to the bifurcation area.In figure 3.3 we compare the numerical tipping to the predicted tipping in (3.17) overa range of . Here we see performance even better than in section section 2.2 as even forrelatively large the prediction has a small error. This is an artifact of having a higherdimensional problem, where now two exponentials in (3.16) are dominating the behavior ofthe solution in the V > 0 region. As in section 2.2, the concavity of the predicted tippingagainst the numerical tipping matches very well and we can expect the prediction to holdfor reasonably small values of .433.2. High Frequency Oscillatory ForcingFigure 3.3: The numerical tipping (red stars) vs. the estimate η2slow (black line) withη1 = 4 and η3 = .375. The top blue line is the tipping of the second eigenvalue. The tippingcriterion on the numerical solution is V > Vsmooth.3.2 High Frequency Oscillatory ForcingWe consider oscillations occurring in the dynamics of the THC that are not originallyencompassed by the Stommel model [1, 9, 14, 17, 18, 22]. We allow η1 to exhibit oscillatorybehavior to account for this. As η1 appears in both equations for V and T , we consideroscillatory forcing on both, but permit their amplitudes to be different. In other words, thecanonical system (3.1) with A,B ∼ O(1), Ω 1 and = 0 which is the purely oscillatoryforcing problem. Under these conditions, like with the one component model in section 2.3,we expect to find oscillations that are attracting. Here stable behavior should act likeoscillations about the equilibria of a reduced inner problem. Thus our analysis must locatethese equilibria as a function of η2 in order to find the bifurcation.To begin our analysis, we note that the dynamics occur on multiple time scales, a slowt and fast R = Ωt. Following a multiple scales approach, we consider V (t) = V (t, R) andT (t) = T (t, R) and substituting this into (3.1), we get the systemVR + Ω−1Vt =Ω−1 (η1 − η2 + η3(T − V )− T − V |V |+A sin(R)) ,TR + Ω−1Tt =Ω−1 (η1 − T (1 + |V |) +B sin(R)) .(3.18)We follow the lower branch to study dynamics near the non-smooth bifurcation. Thus weconsider the system (3.18) with V < 0VR + Ω−1Vt =Ω−1(η1 − η2 + η3(T − V )− T + V 2 +A sin(R)),TR + Ω−1Tt =Ω−1 (η1 − T (1− V ) +B sin(R)) .(3.19)443.2. High Frequency Oscillatory ForcingFrom (3.19), it makes sense to consider an asymptotic expansion for both V and T in termsof the small quantity, Ω−1, withV (t, R) ∼ V0(t, R) + Ω−1V1(t, R) + Ω−2V2(t, R) +O(Ω−3),T (t, R) ∼ T0(t, R) + Ω−1T1(t, R) + Ω−2T2(t, R) +O(Ω−3).(3.20)Substituting (3.20) into the system (3.19) givesV0R + Ω−1V0t + Ω−1V1R + . . . =Ω−1(η1 − η2 + η3(T0 − V0)− T0 + V 20 +A sin(R))+Ω−2(η3(T1 − V1)− T1 + 2V1V0) + . . . ,T0R + Ω−1T0t + Ω−1T1R + . . . =Ω−1(η1 − T0(1− V0) +B sin(R))+Ω−2(−T1(1− V0) + T0V1) + . . . .Here we find the following equations separated by order of Ω−1 withO(1) :{V0R = 0,T0R = 0,(3.21)O(Ω−1) :{V1R + V0t = η1 − η2 + η3(T0 − V0)− T0 + V 20 +A sin(R),T1R + T0t = η1 − T0(1− V0) +B sin(R),(3.22)O(Ω−2) :{V2R + V1t = η3(T1 − V1)− T1 + 2V0V1,T2R + T1t = −T1(1− V0) + T0V1.(3.23)We learn from (3.21) that both our leading order terms are only dependent on the slow time,V0 = V0(t), T0 = T0(t). Much like section 2.3, we must introduce a solvability conditionon the resonant terms to be able to solve for the correction terms, ensuring for the sub-linear growth resulting in a stable solution. Here, we use the Fredholm alternative (2.24)on (3.22)-(3.23) and search for the equilibrium solutions. This is shown in Appendix C.Collecting these solutions leads to the outer solution in original variablesV ∼V0 − Ω−1A cos(Ωt) + . . . ,T ∼T0 − Ω−1B cos(Ωt) + . . . .(3.24)Here both V0 and T0 are the same equilibrium solutions from the static model in theintroduction withT0(V0) =η11− V0 ,0 = η1 − η2 + η3(T0(V0)− V0)− T0(V0) + V 20 .For the one component model in section 2.3, we had to use a local expansion so we neededto scale both the variable x as well as the parameter µ and analyze the local behavior aroundthe axis x = 0. Since we again have non-smooth dynamics at the axis V = 0, we expect touse a local expansion for the two component model as well. While the precise scaling of thebreakdown of the outer solution (3.24) is too complex for us to search for, we instead observethat once V0 → 0 the oscillations begin to dominate the solution and this is not consistent453.2. High Frequency Oscillatory Forcingwith our assumptions of the expansion. So we introduce the local variables analogously tosection 2.3V =Ω−1X,T =η1 + Ω−1Y,η2 =η2ns + Ω−1n.(3.25)Substituting (3.25) into (3.1) leads to the inner systemX˙ =− n+ η3(Y −X)− Y − Ω−1X|X|+ΩA sin(Ωt),Y˙ =− η1|X|−Y − Ω−1|X|Y + ΩA sin(Ωt).(3.26)The form suggests behavior on the same time scales in (3.26), the slow t and the fast R = Ωt.Considering X(t) = X(t, R) and Y (t) = Y (t, R) gives the multiple scales inner systemXR + Ω−1Xt =Ω−1 (−n+ η3(Y −X)− Y )− Ω−2X|X|+A sin(R),YR + Ω−1Yt =Ω−1 (−η1|X|−Y )− Ω−2|X|Y +B sin(R).(3.27)Once more, as we see the small quantity Ω−1 appearing in (3.27), we choose an expansionof the formX(t, R) ∼X0(t, R) + Ω−1X1(t, R) +O(Ω−2),Y (t, R) ∼Y0(t, R) + Ω−1Y1(t, R) +O(Ω−2),(3.28)where we then substitute (3.28) into (3.27) to getX0R + Ω−1X0t + Ω−1X1R + . . . =Ω−1(−n+ η3(Y0 −X0)− Y0) +A sin(R)+Ω−2(X0|X0 + Ω−1X1 + . . . |+η3(Y1 −X1)− Y1) + . . .Y0R + Ω−1Y0t + Ω−1Y1R + . . . =Ω−1(−η1|X0 + Ω−1X1 + . . . |−Y0) +B sin(R)+Ω−2(−|X0 + Ω−1X1 + . . . |−Y0 − Y1) + . . .We then separate the terms by their orders of Ω−1 to find the equationsO(1) :{X0R = A sin(R),Y0R = B sin(R),(3.29)O(Ω−1) :{X1R +X0t = −n− η3X0 − (1− η3)Y0,Y1R + Y0t = −η1|X0|−Y0.(3.30)From (3.29) we find that the leading order terms of (3.28) have the formX0 = P0(t)−A cos(R), Y0 = Q0(t)−B cos(R), (3.31)463.2. High Frequency Oscillatory Forcingwhere Q and P track the slow time components for their respective variables. Substituting(3.31) into (3.30), we apply the Fredholm alternative (2.24) to derive equations for the slowfunctions P0(t) and Q0(t). This givesP0t =− n− η3P0 − (1− η3)Q0,Q0t =−η12pi∫ 2pi0|P0 −A cos(R)|dR−Q0.(3.32)As we are concerned with finding the bifurcation, we search for the equilibrium solutionsto (3.32). We notice that in the equation for Q0 we find a similar integral equation to theinner equation (2.33) in section 2.3. This leads us again to two cases, Case I: |P0(t)|≤ −|A|where the sign of the integrand does not change, and Case II: |P0(t)|< |A| where theintegrand changes and the integral has a non-trivial solution. These cases can be seen infigure 3.4.Figure 3.4: Here we have the parameter ranges for case I and case II shown as the rightmost green vertical line and the bifurcation value as the left blue vertical line respectively.3.2.1 Case I: P0(t) ≤ −|A|We call this the ’below axis’ case, with the solution X0 below the axis V = 0 so we donot expect a bifurcation here. Instead this case helps to simplify the integration in (3.32)but also helps us to determine when the solution acts differently under case II. We use theequilibria of X and Y to define a parameter range in η2 between case I and case II. Underthe conditions of this case, the system (3.32) simplifies toP0t(t) =− n− η3P0(t)− (1− η3)Q0(t),Q0t(t) =η1P0(t)−Q0(t).(3.33)473.2. High Frequency Oscillatory ForcingSolving for the equilibria in (3.33) results inQ0(P0) = η1P0, P0 = − nη1(1− η3) + η3 .Together with these equilibria and with the condition P0(t) ≤ −|A|, we find the param-eter range that distinguishes between case I and case II in terms of our inner parameter,which we then rewrite in original variables withn ≥(η1(1− η3) + η3)|A|,η2 ≥η2ns +(η1(1− η3) + η3)|A|Ω.(3.34)For values of η2 below the values given in (3.34), we see the oscillations crossing above theaxis V = 0 and hence more contribution from the integral in (3.32). We use a separate caseto deal with this behavior.3.2.2 Case II: |P0(t)|< |A|We call this the ’crossing’ case; here the solution oscillates about the axis V = 0 while thecenter of the solution approaches this axis. With this behavior, we expect a bifurcation tooccur in this region and use the equilibria for (3.32) to determine the location. While thisproblem is two-dimensional in nature, the integral in (3.32) is nearly identical to the integral(2.33) in section 2.3. So we may use the ideas of that section here to get an approximatesolution. Thus, under the assumptions of this case, we fix a value of P0 and integrate (3.32)over the regions defined byR1 = arccos(P0/A), R2 = 2pi − arccos(P0/A).As in section 2.3, the solution to (3.32) is negative for P0(t) the region R ∈ [0, R1] andalternates sign for the regions R ∈ (R1, R2] and R ∈ (R2, 2pi]. We follow the same procedureof integrating over each region to get the exact form for (3.32) withP0t =− n− η3P0(t)− (1− η3)Q0,Q0t =−2η1pi(arcsin(P0/A)P0 +√A2 − P 20)−Q0.(3.35)Although this is the explicit inner equation from (3.35), it is analytically too complex tofind an explicit form for a bifurcation and thus we use a second order Taylor approximationto give an approximate systemP0t =− n− η3P0 − (1− η3)Q0,Q0t =−2η1|A|pi−Q0 − η1pi|A|P20 .(3.36)We solve for the equilibria of (3.36) in order to locate the bifurcation. For simplicity, wedefine a = η1pi|A| , and c =2η1|A|pi . Thus the equilibria satisfyQ0(P0) =− aP 20 − c,0 = −n+ (1− η3)c− η3P0 + aP 20 .(3.37)483.2. High Frequency Oscillatory ForcingHere the equation for P0 in (3.37) is a quadratic that would have two solutions, with thesituation of following the lower branch given by the negative solution withP0 =η32a(1− η3) −12a(1− η3)√η23 + 4a(1− η3)(n− c(1− η3)). (3.38)The equilibrium for P0 in (3.38) is real only for positive discriminant. Then the localbifurcation, nosc, is given for vanishing discriminantnosc =η1(1− η3)|A|pi[2−(piη32η1(1− η3))2]. (3.39)Here the equilibria in (3.37) are given in terms of the local variables. We write the solutionfor V , T and bifurcation, η2osc, in the original variablesV (t) ∼Ω−1 (P0 −A cos(Ωt)) ,T (t) ∼η1 − Ω−1(η1pi|A|P20 +2η1|A|pi+B cos(Ωt)),(3.40)η2osc = η2ns +η1(1− η3)|A|piΩ[2−(piη32η1(1− η3))2]. (3.41)With (3.41) we have found the bifurcation induced by the addition of oscillatory forcingin the Stommel model. As we learned from the one component model in section 2.3, theeffect of oscillatory forcing is early bifurcations where η2osc > η2ns. Our result in (3.41)holds under the caveat that we restrict the parameters withη3 <2√2η1pi + 2√2η1,which is the condition to guarantee the second term in (3.41) is positive. This restrictionis reasonable as generally the parameters have the behavior of η3 < 1 and η3 η1 as thethermal variation is much larger in real ocean dynamics than the ratio of relaxation times.In figure 3.5 the numerical solution to (3.1) for V and a zoom of the solution aroundthe numerical bifurcation is shown. The static bifurcation diagram is underlayed as well forcomparison. We contrast the result in (3.41) to these numerics and find that the bifurcationprediction from our analysis agrees. Notice that there is a bifurcation for η2 > η2ns with theoscillations present, this causes a region of the static lower branch in V to never be followed.In figure 3.6 the numerical solution to (3.1) for T and a zoom in around the bifurcation isshown with the static bifurcation diagram underlayed. Due to the bifurcation η2osc > η2nsin V , the maximum value of T is never reached in (b).493.2. High Frequency Oscillatory Forcing(a) (b)(c)Figure 3.5: In (a) the numerical time series solutions to (3.1) is given with parameters ineach qualitatively different case of η2 = {2.3, 1.8, 1.51} with η1 = 4, η3 = .375, A = B = 2and Ω = 10. In (b) these same solutions are shown on the bifurcation diagram. In (c)a zoom in closer to the non-smooth bifurcation region where the blue vertical line is theestimated bifurcation (3.41).503.2. High Frequency Oscillatory Forcing(a) (b)(c)Figure 3.6: In (a) we have the same numerical time series solutions for the qualitativelydifferent cases η2 = {2.3, 1.8, 1.51}. In (b) we plot these solutions over the standard equi-librium plot for V vs. T . In (c) a zoom closer to the bifurcation area is provided.To evaluate the performance of this prediction, we compare (3.41) to the numericalbifurcation over a range of Ω−1. In figure 3.7 we allow for this range to be Ω−1 ∈ (0, .2).For small values, the two agree very well and as we expect, they begin to diverge once thevalues of Ω−1 become too large from the assumption that Ω 1 and the asymptotics cannotcapture the behavior for low frequency oscillations. This outperforms the approximationfrom the one component model in section 2.3.513.2. High Frequency Oscillatory ForcingFigure 3.7: The numerical tipping (red stars) vs. the estimate (black line). The modelparameters are η1 = 4,η3 = .375 and A = B = 2. The bifurcation criterion for thenumerical solution is V > Vsmooth.3.2.3 StabilitySince we have a non-autonomous system when A 6= 0, we approach the stability with alinearized analysis about the equilibria much like in section 2.3. To do this, recall from(3.32) that we found the systemP0t =− n− η3P0 − (1− η3)Q0,Q0t =−η12pi∫ 2pi0|P0 −A cos(R)| dR−Q0.(3.42)We must consider the stability of solutions over the relative sizes of P0(t) with Case I:P0(t) ≤ −|A| and Case II: |P0(t)|≤ |A|.Case I: v0(t) ≤ −|A|For the ’below-axis’ case the solution spends most of its time away from the axis V = 0.We expect to find attraction around the lower branch and thus we expect stability there.Under these conditions, the inner equation (3.42) simplifies to the equationsP0t =−m− η3P0 − (1− η3)Q0,Q0t =η1P0 −Q0.(3.43)The equilibria of (3.43) is found with Q0(P0) = η1P0. Thus we find the following reduced523.2. High Frequency Oscillatory Forcingone component equation and equilibrium Z0 withP0t = −n− (η3 + (1− η3))Q0 = f(P0), Z0 = −nη3 + η1(1− η3) . (3.44)Now we consider a simple linear perturbation about this equilibrium with P0(t) = Z0 +U(t)where ‖U(t)‖ 1. Our standard Taylor expansion about the equilibrium Z0 results inf(P0) =f(Z0) + fP0(Z0)(P0 − Z0) +O((P0 − Z0)2),Ut =− (η3 + η1(1− η3))U(3.45)From (3.45) we now conclude the equilibrium Z0 is hyperbolic and attracting due to theexponential decay in perturbations. Thus we find that no bifurcation occurs for this casewhich agrees with our findings from the analysis above. This holds for η2 in (3.34)η2 ≥ η2ns +(η3 + η1(1− η3))|A|Ω.Case II: |P0(t)|< |A|We called this case the ’crossing’ case and here the solution experiences the non-smoothbehavior when it crosses V = 0. We expect the stability to fail under these conditions andwe found in the analysis that the crossing V = 0 gives (3.42) in the formP0t =− n− η3P0 − (1− η3)Q0,Q0t =−2η1|A|pi− η1pi|A|P20 −Q0.(3.46)As we search for the equilibria of (3.46), we find the equilibrium for Q0 in terms of P0Q0(P0) = −2η1|A|pi− η1pi|A|P20 ,which then gives the following inner equation with the equilibrium Z0 for P0 < 0P0t =− n+2η1|A|pi− η3P0 + η1pi|A|P20 = f(P0),Z0 =pi|A|2η1(1− η3)(η3 −√4η1(1− η3)pi|A| (n− nosc)).(3.47)For simplicity we write Z0 in terms of the local bifurcation nosc we found in the analysiswith (3.39). We now consider a simple linear perturbation about this equilibrium in (3.47)with P0(t) = Z0 + U(t) where ‖U(t)‖ 1. The standard Taylor expansion about theequilibrium is thusf(P0) =f(Z0) + fP0(Z0)(P0 − Z0) +O((P0 − Z0)2),Ut =− 2√η1(1− η3)pi|A| (n− nosc)U.(3.48)533.3. Slow Variation with Oscillatory ForcingThus with (3.48) we learn that the perturbations U(t) decay exponentially as long as thesquare-root is non-zero. We restrict our attention to real solutions so that our perturbationsare real. Thus we have that Z0 is a hyperbolic and asymptotically stable equilibrium forP0 and find stability in Q0 as well. This gives an attracting solution for this region but welose this stability once the square-root becomes zero, here whennosc =η1(1− η3)|A|pi[2−(piη32η1(1− η3))2].This indicates the equilibrium Z0 at this point is non-hyperbolic which is indicative of abifurcation. The results here are in agreement with our analysis and thus we can say thatthe value found in (3.41) is the bifurcation under the oscillatory forcing.3.3 Slow Variation with Oscillatory ForcingWith the one component model solved and both the slowly varying and high oscillatorytwo component models analyzed, we have all of the tools needed to analyze the full system(3.1) with both 1 and A,B ∼ O(1) simultaneously. This is the most general setting wediscuss in this thesis by accounting for both, slowly varying η2 that leads to abrupt changesseen in [1, 14, 17] as well as the oscillatory forcing seen in [19, 9]. Our goal is to study theinteraction of these mechanisms in the physical Stommel model and give an approximationon catastrophic behavior in the model. Under the framework of slowly varying parameterswe expect to find a tipping point instead of a bifurcation. Hence our method for finding thetipping point follows a mixture of both section 3.1 and section 3.2. This procedure dictatesthat we search for inner behavior about the non-smooth bifurcation and to do so we needto solve the inner equation and estimate when this solution abruptly transitions towardsthe upper branch. Ultimately, we provide a solution that captures the abrupt change fromthe lower stable branch to the upper branch in the full Stommel model.To begin, we take our standard approach of following the lower branch towards thenon-smooth behavior with V < 0 in (3.1) which gives the following systemV˙ = η1 − η2 − T + η3(T − V ) + V 2 +A sin(Ωt),T˙ = η1 − T (1− V ) +B sin(Ωt),η˙2 = −.(3.49)As in section 2.4, we write the frequency in terms of the slow variation, Ω = −λ withexponent λ > 0. This assumption allows us to find the relative influence of the slowlyvarying parameter and fast oscillations on the tipping. We notice in (3.49) that there isvariation on a slow scale in η2(t) and on a fast scale in sin(Ωt), so this suggests a multiplescales approach with slow time τ = t and fast time R = −λt. Using this approach in(3.49) yieldsVR + λ+1Vτ = λ(η1 − η2 − T + η3(T − V ) + V 2 +A sin(R)),TR + λ+1Tτ = λ (η1 − T (1− V ) +B sin(R)) ,η2τ = −1.(3.50)543.3. Slow Variation with Oscillatory ForcingTo approach the solution to the outer equations in (3.50), we use an asymptotic expansionwith both λ and integer powers as we have not specified the range of λ and both could besignificant. Thus our expansion isV (τ,R) ∼V0(τ,R) + λV1(τ,R) +O(2λ, λ+1),T (τ,R) ∼T0(τ,R) + λT1(τ,R) +O(2λ, λ+1).(3.51)Here we substitute (3.51) into (3.50) to giveV0R + λ+1V0τ + λV1R + . . . =λ(η1 − η2 − T + η3(T0 − V0) + V 20 +A sin(R))+ 2λ (−η3V1 − (1− η3)T1 + 2V0V1) + . . . ,T0R + λ+1T0τ + λT1R + . . . =λ (η1 − T0(1− V0) +B sin(R))+ 2λ (−T1 + T0V1 + T1V0) + . . . .Separating the equation at each order of then gives the following sets of equationsO(1) :{V0R = 0,T0R = 0,(3.52)O(λ) :{V1R = η1 − η2(τ) + η3(T0 − V0)− T0 + V 20 +A sin(R),T1R = η1 − T0(1− V0) +B sin(R),(3.53)O(2λ) :{V2R + 1−λV0τ = η3(T1 − V1)− T1 + 2V0V1,T2R + 1−λT0τ = −T1(1− V0) + T0V1.(3.54)We consider a λ where the next order in (3.51) is O(2λ). Similar results follow from achoice in λ where O(λ+1) < O(2λ). From (3.52) the leading order terms in our expansionare only dependent on slow time, V0 = V0(τ) and T0 = T0(τ). The solutions of (3.53) and(3.54) are found in Appendix C, giving the outer solutionV ∼V0 + (V0τ (1− V0) + (1− η3)T0τ )(1− η3)T0 + (2V0 − η3)(1− V0) − λA cos(Ωt),T ∼T0 + T0τ1− V0 −T0(V0τ (1− V0) + (1− η3)T0τ )(1− η3)T0(1− V0) + (2V0 − η3)(1− V0)2 − λB cos(Ωt),(3.55)where V0 and T0 are the same leading order solutions from the slowly varying Stommelmodel in section 3.1. Unfortunately, the common theme of the Stommel model is that theseouter solutions are complex but it is clear the outer expansion breaks down as V0 → 0 asthe scale separation between V0 and V1 no longer exists. Thus we derive a scaling for theinner equations which is analogous to that of section 2.4.For simplicity, we assume that the scaling for both V and T is the same, but this isn’tnecessary to arrive at the same conclusion. Hence we rescale about the bifurcation point(η2ns, Vns, Tns) = (η1η3, 0, η1) withV = αX, T = η1 + αY, η2(t) = η2ns + βn(t), (3.56)553.3. Slow Variation with Oscillatory Forcingwhere both α > 0 and β > 0 allow for this to be a local scaling. Applying the local variablesin (3.56) to the full Stommel model (3.1) givesαX˙ =− βn(t)− α(X + (1− η3)Y )− 2αX|X|+A sin(−λt),αY˙ =− α(η1|X|+Y ) + 2α|X|Y +B sin(−λt)n˙ =− 1−β.(3.57)From (3.57) it is apparent that there is still fast behavior within the oscillations. Alsonote that due to the scaling in (3.56), the local variable’s slow behavior has been movedinto the regular time t. To flesh out the particular choice in α, we then take a multiplescales approach to capture the oscillations with t and R = −λt. This choice comes withthe ambiguity in β and is discussed further below. Applying the multiple scales in (3.57)results inα−λXR + αXt =− βn(t)− α(X + (1− η3)Y )− 2αX|X|+A sin(R),α−λYR + αYt =− α(η1|X|+Y )− 2α|X|Y +B sin(R)nt =− 1−β.(3.58)Next we balance the leading order terms in each equation of (3.58), α−λXR and A sin(R)as well as α−λYR with B sin(R), which gives us that α = λ and confirms that the scales Vand T are the same. The scaling for η2 has yet to be determined and could have multiplepossibilities depending on λ. Due to this choice in α we expect the oscillatory term topersist in the inner asymptotic expansion of (3.1) regardless of choice in λ.Here we use a multiple scales approach with t and R = −λt on the full Stommel model(3.1) along with the general scaling (3.56) on η2 which givesVR + λVt =− λ+βn(t)− λ(η1 − η1η3 + η3(T − V )− T − V |V |+A sin(R)),TR + λTt =λ(η1 − T (1 + |V |) +B sin(R)),nt =− 1−β.(3.59)In section 2.4 the results depended on the relative size of the slow variation with respectto the oscillations. The distinction in that case was when λ ≤ 1, where a mixture betweenthe slow variation and oscillations influence the tipping or λ > 1, where the slow variationdominates the tipping. We find a similar distinction here and hence consider a separateasymptotic expansion for Case I: λ ≤ 1 and Case II: λ > 1 to find an accurate classificationof behavior for the full Stommel model.3.3.1 Case I: λ ≤ 1We call this the ’mixed effects’ case where there is significant influence from both slowvariation and fast oscillations due to the size of λ. Here we cannot determine what the nextterm in the expansion should be and thus we choose a general expansion withV (t, R) ∼λX0(t, R) + qX1(t, R) + . . . ,T (t, R) ∼η1 + λY0(t, R) + qY1(t, R) + . . . ,(3.60)563.3. Slow Variation with Oscillatory Forcingwith q > λ to allow for O(λ) to be the leading order term which our local analysis suggested.Substituting (3.60) into (3.59) then gives the governing dynamics for this caseX0R + λX0t + q−λX1R . . . = − βn(t)− λ(η3X0 + (1− η3)Y0)− 2λ(X0 + q−λX1 + . . .)|X0 + q−λX1 + . . . |− q(η3X1 + (1− η3)Y1) +A sin(R) + . . . ,Y0R + λY0t + q−λY1R + . . . =− λ(η1|X0 + q−λX1 + . . . |+Y0 + q−λY1 + . . .)+ 2λ|X0 + q−λX1 + . . . |(Y0 + q−λY1 + . . .)+B sin(R) + . . . .Separating by the distinct orders of gives the following equations at each orderO(1) :{X0R = A sin(R),Y0R = B sin(R),(3.61)O(λ) :{q−2λX1R +X0t = −β−λn(t)− η3X0 − (1− η3)X0,q−2λY1R + Y0t = −η1|X0|−Y0.(3.62)We learn from (3.62) that q = 2λ balances the equations, which implies that λ > 12 foran expansion to be found. If λ ≤ 12 then we would need to include the quadratic terms inO(λ) and our reduced equations would be the same as the full Stommel model (3.1). Thisindicates that our local approximation would no longer hold and that the high frequencyassumption is failing. Without this, the physical behavior of the problem is qualitativelydifferent and we explore this further in chapter 4. On this range of λ, there are two choicesfor the scaling of η2 with β = 1 or β = λ. The advantage to choosing β = λ is that theequation (3.62) is simple, but the slow variation equation has the form nt = −1−λ andimplies a slower time scale. On the other hand, β = 1 keeps a small coefficient on n in(3.62) but gives the slow variation equation as nt = −1, which allows the time scale t to beused. Both of these choices result in the same approximation of the tipping point in originalvariables. Here we choose β = 1 for convenience. From (3.61) we find the appropriate formsof the leading order terms, X0 = P0(t)−A cos(R) and Y0 = Q0(t)−B cos(R). Using theseforms for the leading order term and applying the Fredholm alternative (2.24) to (3.62) wefindP0t =− 1−λn(t)− η3P0 − (1− η3)Q0,Q0t =−η12pi∫ 2pi0|P0(t)−A cos(R)| dR−Q0,nt =− 1.(3.63)Analogously to section 3.2 we must approach this integration with the relative sizeof P0(t) to the amplitude of oscillation A in mind. This is due to the sizes determiningthe contributions of the integral in (3.63). We consider these sizes of P0(t) as Sub-case I:P0(t) ≤ −|A| and Sub-case II: |P0(t)|< |A| similar to section 3.2. These cases separateX0 = P0−A cos(R) from either staying entirely below V = 0 or changing signs respectivelyand allows for the integration to have two distinct forms.573.3. Slow Variation with Oscillatory ForcingSub-Case I: P0(t) ≤ −|A|We call this the ’below-axis’ sub-case, where the solution P0(t) is entirely below the axisV = 0 and this means the full solution X0 has oscillations far from crossing. Under theseconditions we don’t expect any tipping behavior as the solution is far from V = 0, butwe may use this size of P0(t) to find the range of η2 that distinguishes these cases. WithP0(t) ≤ −|A|, we find (3.63) simplifies toP0t =− 1−λn(t)− η3P0 − (1− η3)Q0,Q0t =η1P0 −Q0.(3.64)Although we have the means available to solve (3.64) as it takes the form of an equationwe have seen in section 3.1, instead we search for when the pseudo-equilibrium fails theassumption of this sub-case. This results in the parameter range between these sub-casesand taking this approach is more convenient than solving the system. Here the form of thepseudo-equilibria is simple to find as Q0(P0) = η1P0 and thusP0(t) = −1−λ n(t)η3 + η1(1− η3) .We recall that for this sub-case P0(t) ≤ −|A|, which gives the range for n and we rewritethis in the original variables of η2 withn ≥λ(η3 + η1(1− η3))|A|,η2 ≥η2ns +(η3 + η1(1− η3))|A|Ω.(3.65)With the parameter range (3.65), we now have an effective region for sub-case I and havethe range for sub-case II in terms of the parameter η2.Sub-Case II: |P0(t)|≤ |A|We call this the ’crossing’ sub-case and under these conditions we see the integral in (3.63)is more complex. As this contribution changes, there is an increasing effect on the systemand it is here that we anticipate the tipping point to occur. In section 2.4, we found asimilar integral to (3.63) that we could evaluate with the assumption that A ∼ O(1). Notethat the assumptions that allowed for evaluation of the integral in (2.59) from section 2.4still hold here with a fast time R that is sufficiently large due to the high frequency. Thuswe follow the approach from section 2.4 by integrating with R1 = arccos(P0/A) and R2 =2pi−arccos(P0/A) which is analogous to T1 and T2 from section 2.4 and then use a quadraticTaylor approximation about P0 = 0 as in (2.63) to giveP0t =− 1−λn(t)− η3P0(s)− (1− η3)Q0,Q0t =−2η1|A|pi− η1pi|A|P20 −Q0.(3.66)The form in (3.66) is known as a quadratic two component Riccati-type equation. Werecall the assumption that the solution to the equation for T being in terms of V wasrealistic to the THC, so the behavior we are interested in lies within the dynamics for V .583.3. Slow Variation with Oscillatory ForcingFor this reason, we choose to approximate by reducing (3.66) to a one component model byassuming our equation for Q0 is in pseudo-equilibrium withQ0(P0) = −2η1|A|pi− η1pi|A|P20 . (3.67)The resulting reduced system from introducing the equilibrium (3.67) into the inner equation(3.66) is thenP0t =− 1−λn(t) +2η1(1− η3)|A|pi− η3P0 + η1pi|A|P20 ,nt =− 1.(3.68)Again for convenience, we rewrite the differentiation in (3.68) in terms of the parameter ngivingP0n =1−λn− 2η1(1− η3)|A|pi+ η3P0 − η1(1− η3)pi|A| P20 . (3.69)Now (3.69) is in a form where the result from (1.2) applies. Thus we determine that (3.69) isan Airy-type equation and that the tipping point follows with (1.4). We write the solutionin original variables and obtain results similar to previous sections withntip =− (λ−1)/3(pi|A|η1(1− η3))1/3(2.33810) + λ−1η1(1− η3)|A|pi(2−(piη32η1(1− η3))2),η2tip =(pi|A|η1(1− η3)Ω)1/3µsmooth + η2osc.(3.70)We conclude that the tipping point in (3.70) has a similar form as the tipping pointfound with (2.67) in section 2.4, where we found a weighted average between the smoothtipping point µsmooth from [24] and the oscillatory bifurcation η2osc for this range of λ. Wealso found in (3.62) that any choice for λ ≤ 12 gave equations that cannot be studied usingthe multiple scales approach in this section. This heuristically makes sense since for λ ≤ 12we have low frequency oscillations with our polynomial relationship and the contributions tothe dynamics from this behavior require a different approach than presented in this paper.See [24] for an example of a low-frequency method.3.3.2 Case II: λ > 1We call this case the ’slowly varying dominant’ case. Here we expect integer powers of toappear in leading order due to the O(λ) being quite small for this range of λ and thus wechoose the expansionV (t, R) ∼X0(t, R) + λX1(t, R) + qX2(t, R) + . . . ,T (t, R) ∼Y0(s,R) + λY1(t, R) + qY2(t, R) + . . . .(3.71)Substituting (3.71) into (3.59) then gives593.3. Slow Variation with Oscillatory ForcingX0R + λ+1X0t + λX1R + . . . = − λ+βn(t)− λ+1(η3X0 + (1− η3)Y0)− λ+2(X0 + q−λX1 + . . .)|X0 + q−λX1 + . . . |− 2λ(η3X1 + (1− η3)Y1) + λA sin(R),Y0R + λ+1Y0t + λY1R + . . . =− λ+1(η1|X0 + λ−1X1 + q−1X2 + . . . |−Y0 − λ−1Y1 + . . .)+ 2|X0 + λ−1X1 + q−1X2 + . . . |(Y0 + λ−1Y1 + . . .)+ λB sin(R).We separate by each distinct order of to find the equationsO() :{X0R = 0,Y0R = 0,(3.72)O(λ) :{X1R = A sin(R),Y1R = B sin(R),(3.73)O(λ+1) :{q−λ−1X2R +X0t = −β−1n(t)− η3X0 − (1− η3)Y0,q−λ−1Y2R + Y0t = −η1|λ|X0 + λ−1X1|−Y0.(3.74)We learn in (3.74) that q = λ + 1 balances terms in the expansion. We also find thatβ = 1 gives both simple expressions in (3.74) but also nt = −1. Here a single β is found asopposed to case I where we chose the value of β for convenience. In (3.72) we find that theleading order behavior for this case depends only on slow time, X0 = X0(t) and Y0 = Y0(t),thus giving dominant behavior in this case. With (3.73) we find the correction terms aregiven by X1 = P1(t) − A cos(R) and Y1 = Q1(t) − B cos(R) and since the slow behaviorin X1 and Y1 is just next order corrections to the purely slow X0 and Y0, without loss ofgenerality we set P1 ≡ Q1 ≡ 0. This gives purely oscillatory corrections, X1 = −A cos(R)and Y1 = −B cos(R). Applying Fredholm (2.24) to (3.74) then givesX0t =− n(t)− η3X0 − (1− η3)Y0,Y0t =−η12pi∫ 2pi0|X0(t)− λ−1A cos(R)| dR− Y0,nt =− 1.(3.75)In case I, we used the pseudo-equilibrium of Q0 regardless of the size of the oscillationsto find a solvable equation. Here we expect (3.75) to have a quadratic form like case I,so we choose a priori to use a similar reduction by assuming the equation for Y0 is in it’spseudo-equilibrium withY0(X0) = − η12pi∫ 2pi0|X0 − λ−1A cos(R)| dR. (3.76)603.3. Slow Variation with Oscillatory ForcingWe find the resulting reduced one component equation by introducing (3.76) into the innerequation (3.75) withX0t = −n(t)− η3X0 +η1(1− η3)2pi∫ 2pi0|X0(t)− λ−1A cos(R)| dR. (3.77)The behavior in (3.77) is similar to case I as long as the amplitude of oscillation insidethe integral are consistent with the assumptions of case I, λ−1A ∼ O(1). Under thisassumption, we find that λ ≈ 1 to see mixed behavior of case I and thus we once morefollow the method of section 2.4. Our assumption on the size of the oscillations allow usto integrate (3.77) with R1 = arccos(X0/λ−1A)and R2 = 2pi − arccos(X0/λ−1A)whichagain are analogous to T1 and T2 from section 2.4. Another application of a quadraticTaylor approximation about X0 = 0 as in (2.74) then yieldsX0t = −n(t) + λ−12η1(1− η3)|A|pi− η3X0 + 1−λ η1(1− η3)pi|A| X20 . (3.78)Once more, we find a form to which we apply the result from (1.2). Thus we find thetipping point for the local parameter n with (1.4) and then transform back into the originalvariables to find the tipping point in η2 withnmixed =− (λ−1)/3(pi|A|η1(1− η3))1/3(2.33810) + λ−1η1(1− η3)|A|pi(2−(piη32η1(1− η3))2),η2mixed =(pi|A|η1(1− η3)Ω)1/3µsmooth + η2osc.(3.79)Our result for η2mixed is not surprising as we had similar forms and assumptions to theones of case I. For λ > 1 and away from 1, the amplitude of the oscillations is smaller, andinside the integral in (3.58) the contribution from the oscillations is reduced. Although noexact cut-off exists and, depending on the choice in other model parameters η1 and η3, wefind this is typically for λ ≥ 1.5, the system (3.75) is approximated byX0t =− n(t)− η3X0 − (1− η3)Y0,Y0t =− η3|X0|−Y0,nt =− 1.(3.80)Here (3.80) is the same system as the slowly varying model in section 3.1. With thesame inner equation and slowly varying n, we are able to use the approximation found insection 3.1 for the tipping point η2slow from (3.17) here as well. This indicates that theamplitude of the oscillations are small for larger λ, so that only the slow variation affectsthe tipping of the Stommel model.With both case I and case II, we have described the tipping behavior for any choice in λ.For λ ≤ 1, we found a similar combination of contributions from the oscillatory bifurcationη2osc and the smooth tipping point µsmooth in the Airy equation as in section 2.4. Theweighted averaging is reminiscent of the smooth tipping point analysis in [24] where herethe weight depends on the frequency Ω. With this frequency dependency, we have further613.3. Slow Variation with Oscillatory Forcingevidence to see mixed behavior for medium sized frequency (which corresponds to smallerλ). This behavior is observed for both λ ≤ 1 and λ > 1 but the oscillatory behaviorcontributes less to the advance of the tipping point for larger λ. For λ sufficiently large, theoscillations have a negligible contribution and we recover the tipping point η2slow from theslowly varying model. The results for the tipping point in the Stommel model with slowlyvarying parameter η2 and oscillatory forcing are summarized in the following table.Two Component Tipping Points > 0 and A = 0: η2slow = min(η1η3 − log()/λi) for i ∈ {1, 2} = 0 and A 6= 0 with Ω 1: η2osc = η1η3 + η1(1−η3)|A|piΩ(2−(piη32η1(1−η3))2) > 0, A 6= 0 and 12 < λ ≤ 1: η2mixed =(pi|A|η1(1−η3)Ω)1/3µsmooth + η2osc > 0, A 6= 0 and λ > 1 with λ ≈ 1: η2mixed =(pi|A|η1(1−η3)Ω)1/3µsmooth + η2osc > 0, A 6= 0 and λ > 1: η2slow = min(η1η3 − log()/λi) for i ∈ {1, 2}Table 3.1: Overview of the tipping points in the two component model for each mechanismand case.In figure 3.8, we see an example of the numerical solution of V to the Stommel model(3.1) with slow variation and oscillatory forcing. This example illustrates the tipping forcase I with λ ∈ (12 , 1], so that both the slow variation and oscillatory forcing influence thetipping point. The vertical lines are the tipping, black solid for the numerical and bluedotted for the approximation for this case (3.70). Although there is a mixture of effects, thetipping point gives a value near the oscillatory bifurcation η2osc. This tells us that for thesechoices in the model parameters the strongest effect is the oscillatory forcing. The resultsare shown in figure (3.9) in the V − T plane. Here we see that due to the early tipping inV , the solution for T also never achieves it’s maximum and there is early tipping here aswell, which agrees with the assumptions we had made of considering T responding to V .623.3. Slow Variation with Oscillatory Forcing(a) (b)Figure 3.8: The model values are λ = .8, = .01 with A = B = 2. In (a) the numericalsolution (black dotted line) to (3.1) is given with η1 = 4, η3 = .375. In (b) a zoom in closerto the non-smooth bifurcation region where the blue dotted vertical line is the tippingpoint (3.70) and the black vertical line are the tipping points with the tipping criterionV > Vsmooth on the numerical solution.(a) (b)Figure 3.9: The model values are λ = .8, = .01 with A = B = 2. In (a) we have thenumerical solution (black dotted) over the static equilibrium plot for V vs. T . In (b) azoom of the bifurcation area.In figure 3.10 we have chosen a value of λ in case II as λ > 1 but we also have thatλ ≈ 1. Thus we see comparable behavior to case I with the addition that the slow variationis now dominant. Upon a zoom in, it is apparent that oscillations are still present and wesee a mixture of effects that cause a similar tipping to case I to take place. We’ve plottedthe tipping point η2slow from the slowly varying model (3.17) as the green vertical dottedline for comparison. As the numerical tipping point is moving towards the slowly varying633.3. Slow Variation with Oscillatory Forcingtipping point this confirms that the slow variation is indeed dominating the tipping. Infigure 3.11 we again see very similar behavior to the slowly varying model in section 3.1 butthe zoom-in further reveals the oscillations are present and have minor influence by forcingthe solution to cross the V = 0 axis near the tipping.(a) (b)Figure 3.10: The model values are λ = 1.05, = .01 with A = B = 2. In (a) the numericalsolution (black dotted line) to (3.1) is given with η1 = 4 and η3 = .375. In (b) a zoom incloser to the non-smooth bifurcation region where the blue dotted vertical line is the mixedtipping point (3.70), green dotted verticle line is the slow tipping point (3.17) and the blackvertical line are the tipping points with the tipping criterion V > Vsmooth on the numericalsolution.(a) (b)Figure 3.11: The model values are λ = 1.05, = .01 with A = B = 2. In (a) we have thenumerical solution (black dotted) over the static equilibrium plot for V vs. T . In (b) azoom of the bifurcation area is given.In figure 3.12 we show the numerics for a λ large enough so that the oscillations are643.3. Slow Variation with Oscillatory Forcingnegligible and we recover the slowly varying model in section 3.1. Even upon a zoom it isalmost impossible to see oscillations in this solution. The green dotted vertical line is theslowly varying tipping estimate (3.17) where the blue dotted is the mixed approximation(3.70). Further evidence is seen in figure 3.13, where this figure resembles the slowly varyingV − T plot (3.2).(a) (b)Figure 3.12: The model values are λ = 2, = .01 with A = B = 2. In (a) the numericalsolution (black dotted line) to (3.1) is given with η1 = 4 and η3 = .375. In (b) a zoomin closer to the non-smooth bifurcation region where the blue dotted vertical line is themixed tipping point (3.70), the green dotted vertical line is the slow tipping point (3.17)and the black vertical line are the tipping points with the tipping criterion V > Vsmooth onthe numerical solution.(a) (b)Figure 3.13: The model values are λ = 2, = .01 with A = B = 2. In (a) we have thenumerical solution (black dotted) over the static equilibrium plot for V vs. T . In (b) azoom of the bifurcation area is provided.653.3. Slow Variation with Oscillatory ForcingAlthough the figures above show that we have classified the behavior appropriately forthe various cases in λ and relative solution sizes, performance of this approximate tippingpoint needs to be evaluated to compare with numerical results. In figure 3.14 we comparethe tipping points between case I and case II with the numerically obtained tipping pointsacross a range of λ with a fixed . For smaller λ, the frequency Ω is smaller and theinfluence of the oscillations on tipping become more predominant. Recall the assumptionthat Ω = −λ 1 and that for λ ≤ 12 we observe Ω ∼ O(1). We do not consider lowfrequency corresponding to λ < 12 in this section. For larger λ, there is a reduced influencefor the oscillatory forcing until it is negligible for some λ > 1. We notice that our reductiontipping approximation for λ < 1 has some bias which can be attributed the information lostfrom reducing the full two component Riccati equations to a one component model. Eventhough we use a one component reduced equation to get these approximations, they seemto be performing quite well across all λ which confirms the approach leads to a sufficientapproximation.Figure 3.14: An example of numerical tipping (red stars) as the numerical solution to (3.1)passes the tipping criterion V > Vsmooth. The parameter values are = .01 and A = B = 3.The lines are the case I tipping estimate (black solid line) and the case II tipping estimate(blue dotted line).We also are interested in the performance of the tipping approximations across valuesof for λ fixed, which is seen in figure 3.15. For case I tipping, the range of appropriate is highly dependent on the choice in λ. Often, the range is very small to get accurateestimates which is another artifact of using a reduced model. For case II tipping, we seethat the numerical results are being approximated by the slowly varying tipping section 3.1as the effective oscillations shrink, λ−1A→ 0 as λ→∞.663.3. Slow Variation with Oscillatory Forcing(a) λ = .8 (b) λ = 1.3Figure 3.15: The numerical tipping (red stars) follows the appropriate case depending onλ for = 0.01. The case I tipping estimate η2mixed (black solid line) and slowly varyingtipping estimate η2slow (blue dotted line) are shown.With the numerical results agreeing with our analytic results, we may finally concludethat this method is both useful for analyzing the non-smooth behavior in the Stommelmodel and also results in an approximation that is more accurate in the extremal cases ofthe model (i.e Ω 1 or 1). This gives us a very accessible means of extracting thetipping in the full two component model without needing to solve difficult Riccati equationsor other complex systems that appear in the full problem. We still need to confirm thatthe solutions we’ve found are attracting until stability is lost at the tipping point.3.3.3 StabilityCase I: λ ≤ 1From the analysis, we obtained the inner equations that govern the behavior of the solutionfor this range of λ areP0t =− 1−λn(t)− η3P0 − (1− η3)Q0,Q0t =−η12pi∫ 2pi0|P0 −A cos(R)| dR−Q0.(3.81)But we also found that the relative size of P0(t) dictates the contribution from the integralin (3.81). In the analysis we treat these as Sub-case I: P0(t) ≤ −|A| and Sub-case II:|P0(t)|< |A| that each require a separate analysis.Sub-Case I: P0(t) ≤ −|A|We called this the ’below-axis’ sub-case due to the solution remaining below the axis andthus we anticipate this sub-case to remain attracting to a solution near the lower branch.673.3. Slow Variation with Oscillatory ForcingThe equation (3.81) simplifies for this sub-case toP0t =− 1−λn(t)− η3P0 − (1− η3)Q0,Q0t =η1P0 −Q0.(3.82)From the analysis we choose to reduce (3.82) with the pseudo-equilibria Q0(P0) = η1P0.This gives the following reduced equation with the pseudo-equilibria Z0(t) asP0t =− 1−λn(t)− (η3 + η1(1− η3))P0 = f(t, P0),Z0(t) =− 1−λ n(t)η3 + η1(1− η3) .(3.83)We adopt a similar strategy for analyzing the stability of the one component modelin section 2.4 due to (3.83) being a one-dimensional equation. Hence we perform a simplelinear perturbation about the pseudo-equilibrium with P0(t) = Z0(t)+U(t) and ‖U(t)‖ 1.Taking special care to note that Z0(t) also varies in time, we find the Taylor approximationP0t =f(t, Z0) + fP0(t, Z0)(P0(t)− Z0(t)) +O(‖(P0(t)− Z0(t))2‖2),Ut =− 1−λ 1η3 + η1(1− η3) − (η3 + η1(1− η3))U.(3.84)From (3.84) we find that the perturbations decay exponentially to just below the pseudo-equilibrium Z0. This indicates that the solution for this sub-case is hyperbolically attractingand that there is no tipping for this range of the parameter η2.Sub-Case II: |P0(t)|< |A|We called this the ’crossing’ sub-case and from the analysis above we anticipate the tippingto occur here. The contributions from V > 0 cause the solution to grow and thus weexpect to lose stability. Under the condition |P0(t)|< |A|, we integrate (3.81) with R1 =arccos(P0/A) and R2 = 2pi − arccos(P0/A) and use a Taylor approximation about P0 = 0,which leads toP0t =− n(t)− η3P0 − (1− η3)Q0,Q0t =− λ−12η1|A|pi− 1−λ η1(1− η3)pi|A| P20 −Q0.(3.85)Once more, we assume that Q0 has reached its pseudo-equilibrium to reduce to thefollowing one component inner equation with pseudo-equilibrium Z0(t). Here we let a =η1(1−η3)pi|A| for simplicity, thus the reduced equations areP0t =− 1−λn(t) +2η1(1− η3)|A|pi− η3P0 + aP 20 = f(t, P0),Z0(t) =12a(η3 −√4a(1−λn(t)− nosc)).(3.86)We choose to write the argument of the square root in terms of the local oscillatory bi-furcation nosc found in (3.41). Next, we allow for linear perturbation about the pseudo-equilibrium P0(t) = Z0(t) + U(t) with ‖U(t)‖ 1. Taking a Taylor expansion about the683.3. Slow Variation with Oscillatory Forcingpseudo-equilibrium then allows us to find the local behavior of the perturbations, but re-call that we have contributions in the derivative from the perturbation Ut as well as thepseudo-equilibrium Z0t . This is seen withP0t =Z0t + Ut,Z0t =1−λ√4a(1−λn(t)−nosc)1−λn(t) > nosc,0 1−λn(t) = nosc.(3.87)Thus we find the following Taylor expansion for the perturbationsP0t =f(t, Z0) + fP0(t, Z0)(P0 − Z0) +O(‖P0 − Z0‖2),Ut =1−λ√4a(1−λn(t)−nosc)−(√4a(1−λn(t)− nosc))U 1−λn(t) > nosc,0 1−λn(t) = nosc.(3.88)From (3.88) we find exponentially decaying perturbations to just under the pseudo-equilibriumthis gives asymptotic attraction to a small value for η2 > η2osc which is the oscillatory bi-furcation from section 3.2. This corresponds to loss of the pseudo-equilibrium and so thereis no longer attraction to it, yielding a tipping point for η2 < η2osc which agrees with theresults of our analysis.Case II: λ > 1In the analysis we determined this to be the ’slowly varying dominant’ case and we obtainedthe inner equations that govern the behavior of the solution for this range of λ to beX0t =− n(t)− η3X0 − (1− η3)Y0,Y0t =−η12pi∫ 2pi0|X0 − λ−1A cos(R)| dR− Y0.(3.89)The behavior of this case when λ ≈ 1 is similar to case I, thus we anticipate the similarattraction as well. Hence we consider the behavior when |X0(t)|< λ−1|A| which is thesub-case where we found the tipping point to occur in the analysis. As long as we haveλ−1A ∼ O(1), we are able to follow the same approach as case I where we integrate(3.89) with a R1 = arccos(X0/λ−1A)and R2 = 2pi − arccos(X0/λ−1A)and use a Taylorapproximation about X0 = 0 to getX0t =− n(t)− η2X0 − (1− η3)Y0,Y0t =− λ−12η1(1− η3)|A|pi− 1−λ η1pi|A|X20 − Y0.As in case I, we expect to use the pseudo-equilibrium reduction for Y0 and thus we findthe inner equation with pseudo-equilibrium Z0(t), taking a = η1(1−η3)pi|A| for simplicity,X0t =− n(t) + λ−12η1(1− η3)|A|pi− η3X0 + 1−λaX20 ,Z0(t) =12a(λ−1η3 −√4aλ−1(n(t)− λ−1nosc)).693.3. Slow Variation with Oscillatory ForcingFor a linear stability analysis, we consider the linear perturbation of the pseudo-equilibriumX0(t) = Z0(t) + U(t) with ‖U(t)‖ 1. We apply a Taylor expansion to find a linearequation for perturbations, but again recall that we have contributions to the derivativefrom both the perturbation Ut as well as the pseudo-equilibrium Z0t . This is seen withX0t =Z0t + Ut,Z0t =λ−1√4aλ−1(n(t)−λ−1nosc)n(t) > λ−1nosc,0 n(t) = 1−λnosc.(3.90)Thus we find the following Taylor expansion for the perturbationsX0t = f(t, Z0) + fX0(t, Z0)(X0 − Z0) +O(‖X0 − Z0‖2),Ut =−λ−1√4aλ−1(n(t)−λ−1nosc)−(√λ−14a(n(t)− λ−1nosc))U n(t) > λ−1nosc,0 n(t) = λ−1nosc.(3.91)Similarly to case I, (3.91) shows that the perturbations decay exponentially to justunder the pseudo-equilibrium and we have hyperbolic stability for η2 > η2osc which is theoscillatory bifurcation from section 3.2. After this point is reached, the system loses itspseudo-equilibrium which indicates that the tipping point occurs after the bifurcation η2osc.Comparing this to case I, we see there are small nuances between these perturbations,although the overall stability remains the same. As for when λ gets large, we alreadyestablished that this behaves like the slowly varying model and hence we use the stabilityfrom section 3.1 to conclude that our solution is still stable until the slow tipping pointη2osc.Thus, the stability for both case I and case II agrees with the results found in the analysis.We have that the behavior of our solution is stable from the outer solution and that thisstability holds before the solution begins to cross the axis V = 0. Once the crossing startsto happen, we lose stability at η2osc. Because there is slow variation in this model, there isdelayed behavior and thus the tipping point occurs after the oscillatory bifurcation η2osc. Inboth cases we discovered that the pseudo-equilibrium has a contribution to the derivativeand this in turn causes the perturbations to decay towards a small constant below thepseudo-equilibrium. This means that there is a small region around the pseudo-equilibriumthat attracts the solution.70Chapter 4Summary and Future WorkWith the results found in this thesis, we have accurately described what kind of behavior ispresent about the non-smooth bifurcation when new mechanisms are introduced in both theone component model (2.1) and Stommel model (3.1). The novelty of this work comes fromthe link between delayed tipping methods to the non-smooth Stommel model which thenpaves the way for a more general approach to the broader class of non-smooth dynamicalsystems.To describe a large class of observable behaviors, we considered the mixture of advancedbifurcation due to high oscillatory forcing with frequency Ω 1 and amplitude A ∼ O(1)as well as delayed tipping due to slow variation in the bifurcating parameter at rate 1.We found that addition of these mechanisms have opposite effects on the tipping pointand do mix with a kind of weighted average to produce an effective tipping approximation.The main results of the paper are the relative effects for the non-smooth bifurcation ascompared to the smooth bifurcation. In the case of the one component model, the strengthof the non-smooth influence on the tipping point is vastly different for each component ascompared to the smooth influence. All of the smooth approximations we pull directly from[24] as we delve deeper into a comparison between the two.For the slowly varying parameter with > 0 and no oscillatory forcing with A = 0,we compare the non-smooth tipping point, µslow = ln()/2, to the smooth tipping point,µsmooth = 2/3(−2.3381). Immediately, we notice that these approximations have differentfunctions as seen in figure 4.1. From the figure, it is clear that the delayed non-smootheffects are much smaller and flatten out as compared to the continuously increasing delayedsmooth effects. This indicates that the response to the non-smooth bifurcation occurs soonerthan the response to the smooth bifurcation. Hence we may say that the non-smooth effectsare stronger than the smooth effects.71Chapter 4. Summary and Future WorkFigure 4.1: Comparison between the slow tipping points across . The blue solid line is thenon-smooth tipping points where the red dotted line is the smooth tipping points.For oscillatory forcing with A 6= 0 and static parameter values with = 0, we comparethe non-smooth bifurcation, µosc =4|A|piΩ , to the smooth bifurcation, µsm+osc =A22Ω2. Herewe have a sense of the strength of each case directly from the functions in terms of Ω−1,the non-smooth has a linear response whereas the smooth has a quadratic response. Thisappears clearly in figure 4.2 where we see that the advanced bifurcation in the non-smoothcase is significantly greater than that of the smooth case. This means that the non-smoothbifurcation causes the oscillations to advance the bifurcation much further away as opposedto the smooth bifurcation. Here this also indicates that the effects of the non-smooth effectsare stronger than that of the smooth effects.72Chapter 4. Summary and Future WorkFigure 4.2: Comparison between the oscillatory bifurcation across Ω−1 for A = 2. The bluesolid line is the non-smooth case where the red dotted line is the smooth case.When we combine these mechanisms, we compare the non-smooth tipping point ,µmixed = w(Ω, A)µsmooth + µosc with weight w(Ω, A) = (2|A|piΩ )1/3, to the smooth tippingpoint, µsmooth + µsm+osc. The case where the mixed approximation take the form of µslowfor Ω → ∞ has already been discussed above so we do not consider this here. It is im-portant to note that although we see a similar form between these two, the weight in thenon-smooth tipping point is dependent on both A and Ω so it has significant influence on thevalue for different frequencies. This is shown in figure 4.3 for varying in (a) and Ω varyingin (b). In (a), since the weight w(Ω, A) < 1 for A < O(Ω), then the slope of non-smoothinfluence is smaller than the slope of smooth influence. We also notice the intercept for thenon-smooth influence being significantly larger corresponding to the oscillatory bifurcationµosc being larger than µsm+osc. Together, these indicate that the non-smooth bifurcationhas a stronger advanced tipping while the curves are both positive, and a smaller delayedtipping when both curves are negative. In (b), the effect of the changing weight w(Ω, A) ismost clear and we see the advanced tipping being quite strong for mid-range Ω. It is alsoclear the non-smoothness of the model causes the approximation µmixed to fail for Ω→∞.This case is analogous to when λ > 1 from the one component analysis and allows us torefer back to the slowly varying problem with > 0 and A = 0 for Ω→∞. For both plots,we conclude the non-smooth effects are stronger than the smooth effects even under themixed case.73Chapter 4. Summary and Future Work(a) (b)Figure 4.3: Comparison between the mixed tipping in (a) with a fixed Ω−1 = .1 and (b)with a fixed = .1. The blue solid line is the non-smooth case where the red dotted line isthe smooth case.These results give insight into the hysteresis behavior of the Stommel model and the lessstudied realm of non-smooth dynamics. The main approach used asymptotic expansions aswell as the methods of multiple scales to identify reduced equations and to find asymptoticsolutions to the models. We found that with oscillatory forcing, the reduced equationshave an expression dependent on the relative size of the solution to the amplitude of oscil-lation. In the smooth version of this problem, this type of behavior was not present andconsidering a case-by-case argument was not necessary. We also discovered that linking theslow variation and the frequency Ω gives important insight into how the system behaves.The methods developed for finding tipping points in the one component model (2.1) gavegood approximations with the numerical results. Due to the many similarities to the twocomponent system (3.1) we were able to modify the same analysis to find the tipping pointshere as well. We also anticipate the non-smooth bifurcation of the Stommel model to havea stronger influence on the solution than the smooth bifurcation.Future work would need to be done on cases where Ω ∼ O(1) or smaller. This case isqualitatively different as slow oscillations may have dramatic contributions to the dynamics.These effects are also seen from the analysis where low frequency oscillations no longer allowfor asymptotic expansions in terms of Ω−1 and no longer fall under our assumptions tointegrate with T1 and T2. Hence this case can influence tipping in a way we hadn’t exploredin this paper, we show an example of this in figure 4.4. Also, large amplitude behaviorA > O(1) can force an additional rescaling before any familiar approaches hold which isseen in figure 4.5. These cases were mentioned but have yet to be performed on this model,although both have been studied around the smooth case in [24]. It is possible that theycould have some surprising results in the non-smooth case. Together, these could help tofurther classify the tipping behavior for the variety of cases in real world ocean dynamics.74Chapter 4. Summary and Future Work(a) (b)Figure 4.4: Low Frequency: Model parameters are = .01, A = B = 1 and Ω = 3.(a) (b)Figure 4.5: Large Amplitude: Model parameters are = .01, A = B = 300 and Ω = 1000.Lastly, we considered only deterministic behavior throughout this analysis but there aremany reasons to incorporate stochastic elements into the Stommel model as well, see [13].From [24] it is concluded that stochastic forcing has elements of both early bifurcationsand delayed tipping and thus a natural follow-up to the analysis in this thesis. We couldconsider stochastic forcing withV˙ = η1 − η2 + η3(T − V )− T − V |V |+Aξ1(t),T˙ = η1 − T (1 + |V |) +Bξ2(t),η˙2 = −V (0) =V 0, T (0) = T 0, η2(0) = η20,(4.1)where ξi(t) is standard Gaussian noise with mean 0 and variance t and initial conditions75Chapter 4. Summary and Future Workcentered on the lower branch. This is shown in figure 4.6 and it is clear that a completelyseparate analysis is needed.(a) (b)(c) (d)Figure 4.6: Stochastic: In (a) many realizations of the numerical solution for V in (4.1) aregiven with model parameters η1 = 4, η3 = .375, = .01 and A = B = .7. In (b) a zoomin closer to the non-smooth bifurcation region can be seen. In (c) we have the realizationsover the standard equilibrium plot for V vs. T . In (d) a zoom of the bifurcation area isshown.76Bibliography[1] Richard B Alley, Jochem Marotzke, William D Nordhaus, Jonathan T Overpeck,Dorothy M Peteet, Roger A Pielke, RT Pierrehumbert, PB Rhines, TF Stocker, LD Tal-ley, et al. Abrupt climate change. science, 299(5615):2005–2010, 2003.[2] David Angeli, James E Ferrell, and Eduardo D Sontag. Detection of multistability,bifurcations, and hysteresis in a large class of biological positive-feedback systems.Proceedings of the National Academy of Sciences, 101(7):1822–1827, 2004.[3] Alain Bensoussan, Jacques-Louis Lions, and George Papanicolaou. Asymptotic analysisfor periodic structures, volume 374. American Mathematical Soc., 2011.[4] Levke Caesar, Stefan Rahmstorf, Alexander Robinson, Georg Feulner, and V Saba.Observed fingerprint of a weakening atlantic ocean overturning circulation. Nature,556(1):191–196, 2018.[5] Henk A Dijkstra. Nonlinear climate dynamics. Cambridge University Press, 2013.[6] T Erneux and JP Laplante. Jump transition due to a time-dependent bifurcationparameter in the bistable ioadate–arsenous acid reaction. The Journal of ChemicalPhysics, 90(11):6129–6134, 1989.[7] Richard Haberman. Slowly varying jump and transition phenomena associated withalgebraic bifurcation problems. SIAM Journal on Applied Mathematics, 37(1):69–106,1979.[8] Angela Hohl, HJC Van der Linden, Rajarshi Roy, Guillermo Goldsztein, FernandoBroner, and Steven H Strogatz. Scaling laws for dynamical hysteresis in a multidimen-sional laser system. Physical review letters, 74(12):2220, 1995.[9] Peter Huybers and Carl Wunsch. Obliquity pacing of the late pleistocene glacial ter-minations. Nature, 434(7032):491, 2005.[10] Amitabh Joshi, Wenge Yang, and Min Xiao. Dynamical hysteresis in a three-levelatomic system. Optics letters, 30(8):905–907, 2005.[11] Peter Jung, George Gray, Rajarshi Roy, and Paul Mandel. Scaling law for dynamicalhysteresis. Physical review letters, 65(15):1873, 1990.[12] Yuri A Kuznetsov. Saddle-node bifurcation. Scholarpedia, 1(10):1859, 2006.[13] MN Lorenzo, JJ Taboada, and I Iglesias. The role of stochastic forcing in climatemodels: The case of thermohaline circulation. In Climate Models. InTech, 2012.77[14] Jochem Marotzke. Abrupt climate change and thermohaline circulation: Mechanismsand predictability. Proceedings of the National Academy of Sciences, 97(4):1347–1350,2000.[15] James Dickson Murray. Asymptotic analysis, volume 48. Springer Science & BusinessMedia, 2012.[16] W. Park and M. Latif. Atlantic meridional overturning circulation response to idealizedexternal forcing. Climate Dynamics, 39(7-8):1709–1726, 10 2012.[17] Stefan Rahmstorf. The thermohaline ocean circulation: A system with dangerousthresholds? Climatic Change, 46(3):247–256, 2000.[18] Stefan Rahmstorf. Ocean circulation and climate during the past 120,000 years. Nature,419(6903):207–214, 2002.[19] Andrew Roberts and Raj Saha. Relaxation oscillations in an idealized ocean circulationmodel. Climate Dynamics, 48(7-8):2123–2134, 2017.[20] Nestor E Sanchez. The method of multiple scales: asymptotic solutions and normalforms for nonlinear oscillatory problems. Journal of Symbolic Computation, 21(2):245–252, 1996.[21] Ru¨diger Seydel. Practical bifurcation and stability analysis, volume 5. Springer Science& Business Media, 2009.[22] M Stastna and WR Peltier. On box models of the north atlantic thermohaline circula-tion: Intrinsic and extrinsic millennial timescale variability in response to deterministicand stochastic forcing. Journal of Geophysical Research: Oceans, 112(C10), 2007.[23] Henry Stommel. Thermohaline convection with two stable regimes of flow. Tellus,13(2):224–230, 1961.[24] Jielin Zhu, Rachel Kuske, and Thomas Erneux. Tipping points near a delayed saddlenode bifurcation with periodic forcing. SIAM journal on applied dynamical systems,14(4):2030–2068, 2015.78Appendix AThe Stommel ModelThe original Stommel model has the following balance equations [23]B1dT1dt=CT1 (Ts1 − T1) + |V |(T2 − T1),B2dT2dt=CT2 (Ts2 − T2) + |V |(T1 − T2),B1dS1dt=CS1 (Ss1 − S1) + |V |(S2 − S1),B2dS2dt=CS2 (Ss2 − S2) + |V |(S1 − T2),(A.1)where B1 and B2 are the volumes of each box containing water of temperature and salinity(T1, S1) and (T2, S2) respectively. The surface conditions of each box are then Tsi and Ssiand the relaxation coefficients for the temperature and salinity are CT and CS . V is theflow rate between the boxes and is linearly related to the density difference. To reduce thisto a two dimensional model that contains all the same information, we use the differencesT = T1 − T2 and S = S1 − S2. If we also assume the relaxation coefficient is proportionalto volume, then we have CTB1= CTB2= RT andCSB1= CSB2= RS . This then gives (A.1) asdTdt=RT (Ts1 − T s2 − T )− |V |T(1B1− 1B2),dSdt=RS(Ss1 − Ss2 − S)− |V |TB1 +B2B1B2.(A.2)To simplify further we scale time, temperature, salinity and the flow rate byt← 1CTt,T ← B1B2CTγαT (B1 +B2)T,S ← B1B2CTγαS(B1 +B2)S,V ←B1B2RTB1 +B2V,(A.3)where αT and αS are the thermal expansion and saline contraction coefficients and γ is a79Appendix A. The Stommel Modelhydraulic constant. Then using (A.3) in (A.2) givesdTdt=(T s1 − T s2 )γαT (B1 +B2)B1B2RT− T − |V |T,dSdt=(Ss1 − Ss2)γαS(B1 +B2)RSB1B2R2T− RSRTS − |V |S.(A.4)Here we define the non-dimensionalized parameters• η1 - the strength of thermal forcing,η1 =(T s1 − T s2 )γαT (B1 +B2)B1B2RT.• η2 - the strength of freshwater forcing,η2 =(Ss1 − Ss2)γαS(B1 +B2)RSB1B2R2T.• η3 - the ratio of the freshwater to thermal relaxation time scales,η3 =RSRT.We recall that V was related to the difference in temperature and salinity, thus givingV = T − S. This results in the non-dimensionalized Stommel modeldTdt=η1 − T − |T − S|T,dSdt=η2 − η3S − |T − S|S.(A.5)It is more convenient to analyze ocean circulation in terms of the flow rate V which is alsoknown as the circulation strength. Then using V = T − S in (A.5) givesdTdt=η1 − T (1 + |V |),dVdt=η1 − η2 − η3(T − V )− T − V |V |.(A.6)80Appendix BOne ComponentHigh Frequency Oscillatory ForcingHere we continue the analysis to explicitly find the solution of the outer equation for thepurely oscillatory model. Recall that we found x1 = v1(t) − A cos(T ), we then apply theFredholm alternative (2.24) to the O(Ω−2) equation in (2.23) to get0 =12pi∫ 2pi0−x1t − 2x1 + 2x0x1 dT,v1t =− 2v1 + 2(1−√1 + µ)v1,v1t =− 2√1 + µv1.(B.1)We search for the equilibrium to find stable behavior on this order but since (B.1) has avery simple form, the equilibrium is v1(t) ≡ 0 and thus we find the correction term to onlyhave oscillatory behavior, x1 = −A cos(T ).Slow Variation and Oscillatory ForcingHere we continue to find the terms of the outer solution for the slowly varying and oscil-latory forcing model. We have thus far found x0 = x0(τ) and we have equations at O(λ)and O(2λ) that give information about x0 and x1 respectively. We apply the Fredholmalternative (2.24) to the O(λ) equation (2.49) to find0 =12pi∫ 2pi0−µ(τ)− 2x0(τ) + x0(τ)2 +A sin(T ) dT,0 =− µ(τ)− 2x0(τ) + x0(τ)2,x0(τ) =1−√1 + µ(τ),x1T =A sin(T ).(B.2)From (B.2) we find that x1 = v1(τ) − A cos(T ), which gives us access to solving the nextorder equation. Thus we now do the same for the O(2λ) equation (2.50) to find0 =12pi∫ 2pi0−1−λx0τ − 2x1 + 2x0x1 dT,1−λx0τ =− 2v1 + 2(1−√1 + µ(τ))v1,v1(τ) =− 1−λ x0τ2√1 + µ(τ).(B.3)81Slow Variation and Oscillatory ForcingRecall that µτ = −1 and that x0τ = − µτ2√1+µ(τ)= 12√1+µ(τ)and thus we find the form ofthe next order term in the expansion asx1(τ, T ) = −1−λ 14(1 + µ(τ))−A cos(T ). (B.4)82Appendix CTwo ComponentHigh Frequency Oscillatory ForcingHere we show that the correction term of the outer solution is purely oscillatory. Fromthe analysis, we found both leading order terms to be purely slow time dependent, i.e.V0 = V0(τ) and T0 = T0(τ). To find the explicit form for these, we apply Fredholm (2.24)to the O(Ω−1) equations in (3.22) to find{0 = 12pi∫ 2pi0 −V0t + η1 − η2 + η3(T0 − V0)− T0 + V 20 +A sin(R) dR,0 = 12pi∫ 2pi0 −T0t + η1 − T0(1− V0) +B sin(R) dR,{V0t = η1 − η2 + η3(T0 − V0)− T0 + V 20 ,T0t = η1 − T0(1− V0),V1R = A sin(R), T1R = B cos(R).(C.1)Since we have a fixed parameter η2, we find the equilibria V0 and T0 as well as the form ofthe correction termsT0(V0) =η11− V0 ,0 =η1 − η2 + η3(T0(V0)− V0)− T0(V0) + V 20 ,V1 =X1(t)−A cos(R), T1 = Y1(t)−B cos(R).Here we note these equilibria to be the same as in the static problem with no forcing fromthe introduction. But with the form of the correction terms, we now solve the equation atO(Ω−2) (3.23) by applying Fredholm (2.24). This results in{0 = 12pi∫ 2pi0 −V1t + η3(T1 − V1)− T1 + 2V0V1 dR,0 = 12pi∫ 2pi0 −T1t + T1(1− V0) + T0V1 dR,{X1t = η3(Y1 −X1)− Y1 + 2X0X1,Y1t = Y1(1−X0) + Y0X1.(C.2)We then search for the equilibria of (C.2) and findY1(X1) = − Y0X11−X0 ,0 =(η3(Y01−X0 − 1)− Y01−X0 + 2X0)X1.Thus we find that the correction terms are purely oscillatory since X1 ≡ 0 and Y1 ≡ 0. Thisgives V1 = −A cos(R) and T1 = −B cos(R).83Slow Variation and Oscillatory ForcingSlow Variation and Oscillatory ForcingHere we continue to find terms of the outer solution by working through the equations(3.53)-(3.54). In the analysis, we had already determined that the leading order terms arepurely slow time dependent, V0 = V0(τ) and T0 = T0(τ). To find their exact form, we applyFredholm (2.24) to the O(λ) equation (3.53) to find{0 = 12pi∫ 2pi0 η1 − η2(τ) + η3(T0 − V0)− T0 + V 20 +A sin(R) dR,0 = 12pi∫ 2pi0 η1 − T0(1− V0) +B sin(R) dR,{0 = η1 − η2(τ) + η3(T0 − V0)− T0 + V 20 ,0 = η1 − T0(1− V0),V1R = A sin(R), T1R = B cos(R).(C.3)The leading order solution to (C.3) is the same as the slowly varying problem from section 3.1withT0(V0) =η31− V0 ,0 =η1 − η2(τ) + η3(T0(V0)− V0)− T0(V0) + V 20 .We also find the form of the correction terms, V1 = X1(τ) − A cos(R) and T1 = Y1(τ) −B cos(R), which allow us to solve the O(2λ) equation (3.54). Applying Fredholm (2.24)here results in{0 = 12pi∫ 2pi0(−1−λV0τ + η3(T1 − V1)− T1 + 2V0V1 +A sin(R)) dR,0 = 12pi∫ 2pi0(1−λT0τ − T1(1− V0) + T0V1)dR,{1−λV0τ = η3(Y1 −X1)− Y1 + 2V0X1,1−λT0τ = Y1(1− V0) + T0X1.(C.4)Recalling that η2τ = −1 and solving (C.4) requires the derivatives of V0 and T0 which aresolvable explicitly asT0τ (V0τ ) = −η1V0τ1− V0 ,V0τ =(1− V0)η1 + η3(η1 + 1− V0)− 2V0(1− V0) .With everything put together, we find the solution to (C.4) asY1(X1) =1−λT0τ − T0X11− V0 ,X1 =1−λ(V0τ (1− V0) + (1− η3)T0τ )(1− η3)T0 + (2V0 − η3)(1− V0) .Here we now have the first correction term as84Slow Variation and Oscillatory ForcingV1(τ,R) =1−λ(V0τ (1− V0) + (1− η3)T0τ )(1− η3)T0 + ((2V0 − η3)(1− V0) −A cos(R),T1(τ,R) =1−λT0τ1− V0 −1−λT0(V0τ (1− V0) + (1− η3)T0τ )(1− η3)T0(1− V0) + (2V0 − η3)(1− V0)2 −B cos(R).85
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Non-smooth dynamics in the Stommel model of the thermohaline...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Non-smooth dynamics in the Stommel model of the thermohaline current Griffith, Cody 2018
pdf
Page Metadata
Item Metadata
Title | Non-smooth dynamics in the Stommel model of the thermohaline current |
Creator |
Griffith, Cody |
Publisher | University of British Columbia |
Date Issued | 2018 |
Description | We analyzed the non-smooth dynamics for the Stommel model for thermohaline circulation with additional mechanisms like slowly varying bifurcation parameters and high frequency oscillatory forcing. Our goal was to find an analytic approximation to the tipping point and forced bifurcation induced by the new features in the model. We first analyze a simpler one component model that has similar structure to the Stommel model and gradually build in more complexity into the one component problem and study the effects on the tipping point.In this context, we compare the relative strengths of the non-smooth effects on the tipping point to the smooth effects which has been previously studied. With the one component model understood, we then apply similar methods to the Stommel model and study thee ffects on the non-smooth tipping points and forced bifurcations. With these results we have the ability to fully describe the hysteresis found in the Stommel model. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2018-08-20 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0371123 |
URI | http://hdl.handle.net/2429/66825 |
Degree |
Master of Science - MSc |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2018-09 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2018_september_griffith_cody.pdf [ 3.01MB ]
- Metadata
- JSON: 24-1.0371123.json
- JSON-LD: 24-1.0371123-ld.json
- RDF/XML (Pretty): 24-1.0371123-rdf.xml
- RDF/JSON: 24-1.0371123-rdf.json
- Turtle: 24-1.0371123-turtle.txt
- N-Triples: 24-1.0371123-rdf-ntriples.txt
- Original Record: 24-1.0371123-source.json
- Full Text
- 24-1.0371123-fulltext.txt
- Citation
- 24-1.0371123.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0371123/manifest