Study of Poincaré Map and limit cycles for Non-smooth Welander’s Ocean Convection Model

Ocean circulation plays a very important role in Earth’s climate system. For instance, heat and moisture are distributed around the globe through ocean circulation. Even small changes, in its patterns, can have considerable effects in both regional and global climates. Because of this, understanding this phenomenon is crucial in predicting events that may become catastrophic. In this sense, we can cite the modern phenomenon “El Niño” and the paleo event “Dansgaard Oeschger”. These events can be linked to changing ocean dynamics, causing global changes in temperature, see [9, 10, 29]. For example, after the Last Glacial Maximum (LGM), also referred the Last Glacial Coldest Period, it has been observed Dansgaard Oeschger events, which are periodic weakening and strengthening of the Atlantic Meridional Overturning Circulation (AMOC). It is believed that those oscillations are caused by a sudden influx of freshwater into the Atlantic Ocean, and such a phenomenon can be well modeled as a relaxation oscillation in a low-dimensional model of temperature and salinity in the AMOC. On the other hand, when paying attention to the Gulf Stream, in a shorter time scale, there is a periodic variability of smaller magnitude in its strength, see [6, 20, 31]. Since the Gulf Stream is a major transporter of heat to Europe, changes in this ocean circulation may imply consequences for systems relying on regular climate patterns, such as agriculture. Therefore, the importance of the role played by ocean circulation in the dynamics of the global climate is undeniable. Then to model and explain any of its variations is obviously extremely important to the world, and for example, the book by Kaper and Engle [21] gives a good introduction to mathematical approaches in climate modeling.

In 1982, Welander [32] introduced a simple box model illustrating that interactions in salinity and temperature between the shallow and deep ocean can produce self-oscillations even with no external forcing, justifying its importance in the climate community. We are particularly interested in Welander’s model since it is an example of a Fillipov system [11] wherein the discontinuity set is created by a varying density causing the ocean to oscillate between stratification and mixing. Such a model is an ocean circulation box model, splitting the ocean into boxes, see Figure 1, and this kind of model has been used to study small and large-scale changes in AMOC. Assuming that the water in each box is well mixed, ordinary differential equations describe the flow of water from one box to the other, and the only relevant variables are the temperature and salinity of the water in each box. In addition, the density is described as a linear combination of the variables.

Observe that the model of Welander separates the ocean into two regions surface and deep ocean. The interaction between those regions is dictated by the relative density of the water, controlled by the temperature and salinity through a linear equation of state $\rho=-\alpha_{1}T+\gamma S$. We assume that the deep ocean have constant density $T_{0}=S_{0}=\rho_{0}=0$. The surface is considered well mixed, where the circulation given by the equations must represent a large-scale behavior such that the density does not have an important impact. Then the equations for the temperature and salinity are given by

where $T_{A}$ and $S_{A}$ are generic atmospheric forcing caused by precipitation, evaporation, and solar forcing, which are considered constants. In addition, $\alpha_{1},\gamma,k_{T}$ and $k_{S}$ are also assumed constant.

The main feature of Welander’s model is the convective mixing nonnegative function $\overline{k}(\rho)$, which could be smooth or piecewise smooth. As the ocean is in general highly stratified, without much mixing between the layers, Welander in his model assumes that convection is large when the density of the surface ocean is large, small when the density is small, and the transition between these states occurs abruptly. These aspects are feasible assumptions, given that the ocean is highly stratified, not having much mixing between its layers, and since the exchange between the layers is considered to be fully turbulent, the same $k$ function is used for temperature and salinity. Therefore, the switching behavior allows some periodicity to appear in the model. To obtain a more simple system, we make a change of coordinates and time rescale given by

such that, renaming the parameters, it is possible to obtain the the following planar differential system

where “$\;\dot{}\;$” denotes the derivative in relation to the time, and the new parameters are $\beta=k_{S}/k_{T}$, $\alpha=(T_{A}\alpha_{1})/(S_{A}\gamma)$ and $k(\rho)=\overline{k}(\rho)/k_{T}.$

For some parameter values, Welander showed an oscillatory movement in its system, using two specific functions $k(\rho)$. First, he chooses a continuous function described by

and numerically he finds a stable periodic orbit. In this work, we are interested in studying Welander’s model [32, Sec. 4] for the other function used by him, when $k(\rho)$ is piecewise smooth given by

where $\varepsilon$ consist of a small parameter, $k_{0}$ is a small constant or zero and, $k_{1}$ is a large constant. Using the function (4), Welander does not discuss the mathematics behind this model with rigor, but also numerically finds a periodic solution for specific values of parameters. However, the mathematical mechanism behind this oscillation is different in the non-smooth model.

As we will see in the course of the work the oscillations in the non-smooth case are caused by the fact that, for some parameter values, stable equilibrium points do not exist in the phase portrait (called virtual equilibrium points). In this way, when the system approaches such a virtual equilibrium, it crosses the discontinuity line and reverses its direction to approach the other virtual equilibrium point. We will see that this behavior can lead to a stable routine.

Therefore, our main objective here is to study the occurrences or not of these movements in an analytical way, for different values of parameters. More specifically, we would like to study the existence of isolated periodic orbits, so-called limit cycles. To keep the paper to a reasonable length, we restrict our attention to the case in which $\beta>0$.

For planar piecewise smooth differential systems we have two kinds of limit cycles: sliding limit cycles and crossing limit cycles. The first contains some segment of the discontinuity set, and the second one only contains some points of the discontinuity set. For more details on these kinds of differential systems see the books [1, 11, 26, 30]. Then, here, we are going to study crossing limit cycles, only mentioned as limit cycles, surrounding some sliding segment of the discontinuity line, for the piecewise smooth Welander’s model (2) considering (4).

Therefore, the piecewise smooth Welander’s system can be written as a linear planar piecewise smooth system,

where the discontinuity region (or discontinuity set) $\Sigma$ is given by $f^{-1}(0)$, with $0$ being a regular value of $f$.

The study of a maximum quota of limit cycles for systems of kind (5) has been the object of study for many years. Lum and Chuca in [25], under the continuity hypothesis conjectured that system (5) had at most one limit cycle. In 1998, Freire et al in [12] give the first proof for such conjecture. In 2021, Carmona et al [5], using new techniques, provide a new more simple proof for this conjecture. Han and Zhang in [17] gave the first example of a piecewise linear differential system with two limit cycles. And, in the year of 2012, Huan and Yang [18] using a numerical argument showed a system with three limit cycles. However, LLibre and Ponce [24] showed a first analytical proof of a linear piecewise smooth system with three limit cycles. Actually, there are many works that provided examples with three limit cycles for this kind of differential system, see [3, 4, 14, 15, 27, 28].

Using the theory developed in [19], remember that we are assuming $\beta>0$, and defining the following new parameters

our main results are:

This paper is organized as follows. In Section 2, we present preliminary concepts and results, and we show a normal form for Weelander’s system. The Section 3 is dedicated to developing the Poincaré Map and studying its properties. Finally, in Section 4, we prove Theorems 1 and 2, and we exhibit a non-smooth example of the Welander model with a stable limit cycle, which is the point limit of the smooth example computed by Welander with the function (3).

In this section, we review some definitions and make some accurate analyses of the piecewise smooth Welander’s system that will be used to prove our main results. First, we will review some definitions about piecewise smooth systems. A piecewise smooth system is a system given by

where $f:\mathbb{R}^{2}\rightarrow\mathbb{R}$ is a $C^{1}$ function such that $0$ is a regular value. The discontinuity curve is given by $\Sigma=f^{-1}(0)$, and $Z^{\pm}=(X^{\pm},Y^{\pm})$. Then, the manifold $\Sigma$ divides the plane into two half-planes, $\Sigma^{\pm}$ and the trajectories on $\Sigma$ are defined as follows.

Given a point $p\in\Sigma$, we say that $p$ is a crossing point if, and only if $Z^{+}f(p)\cdot Z^{-}f(p)>0$ where $Z^{\pm}f(p)=\langle\nabla f(p),Z^{\pm}(p)\rangle$ is the Lie Derivative. Moreover, we have $p$ a positive (negative) crossing when $Z^{-}f(p)>0$ and $Z^{+}f(p)>0$ ($Z^{-}f(p)<0$ and $Z^{+}f(p)<0$). And we say that $p$ is sliding (escaping) point if $Z^{+}f(p)\cdot Z^{-}f(p)<0$ such that $Z^{+}f(p)<0$ and $Z^{-}f(p)>0$ ($Z^{+}f(p)>0$ and $Z^{-}f(p)<0$). For crossing points, the solutions in relation to $\Sigma$, approach on one side and leave on the other, and so, its solutions are given by the concatenation of the solutions in $\Sigma^{\pm}$. On the other hand, for sliding or escaping points, the solutions approach or push away of both sides of $\Sigma$, then we need a solution which is continued for future or past in the discontinuity set, making sense. In this way, Fillipov proposed sliding dynamics for this case, which is governed by the Filippov sliding vector field,

and $\lambda$ is such that $Z^{s}f(p)=0$. Therefore, either sliding or escaping points follow the solutions of the Filippov sliding vector field in the appropriate direction. For more details see [11] and, the configuration of the points $p\in\Sigma$, mentioned before, are illustrated in Figure 2.

It is noteworthy that there are other different types of points in the discontinuity set $\Sigma$. For example, we say that $p\in\Sigma$ is a fold point of order $n$ in $Z^{+}$ ($Z^{-}$) with $n\in\mathbb{N}$, $n$ even, if $Z^{+}f(p)=0$ and $(Z^{+})^{n}f(p)\neq 0$ ($Z^{-}f(p)=0$ and $(Z^{-})^{n}f(p)\neq 0$), where $(Z^{\pm})^{n}f(p)=\langle\nabla(Z^{\pm})^{n-1}(p),Z^{\pm}(p)\rangle$ for $n\geq 2$. The fold point for the vector field $Z^{+}$ ($Z^{-}$), is said visible or invisible, respectively, if $(Z^{+})^{n}f(p)>0\,((Z^{-})^{n}f(p)<0)\;\mbox{or}\;(Z^{+})^{n}f(p)<0\,((Z^{-})^{n}f(p)>0).$

Moreover, we say that $p\in\Sigma$ is a focus-focus type if both systems $Z^{+}$ and $Z^{-}$ have an equilibrium point at $p$ with complex eigenvalues. A point $p\in\Sigma$ is a focus-parabolic type if: the system defined in $Z^{+}$ has an equilibrium point of focus type at $p$ while the system $Z^{-}$ has a fold point at $p$, the case parabolic-focus is analogous. Finally, we say that $p\in\Sigma$ is a parabolic-parabolic type if the systems $Z^{+}$ and $Z^{-}$ have a fold point at $p$.

Clearly, the piecewise-smooth system (6) inherits the equilibria of $Z^{+}$ and $Z^{-}$. However, the equilibrium points may or may not appear in the phase portrait, depending on their location in the plane. When they appear, we say they are real, when they don’t, we say they are virtual. We can also pay some attention to other equilibrium points of equilibria of the system (6). It is not difficult to see that $p$, being a sliding or escaping point with $Z^{+}(p)$ and $Z^{-}$ vectors dependent linearly, vanishes $Z^{s}(p)$ in (2). These kinds of points are pseudo-equilibrium points and are internal points of a sliding segment, see [11, 16]. And, an equilibrium of the vector fields $Z^{\pm}$ belonging to $\Sigma$, is called a boundary equilibrium.

In what follows in this section, we will turn our attention to analyzing the non-smooth version of Welander’s system (2) considering the non smooth function (4).

Then, observe that the splitting manifold is given by $\{\rho=-\alpha T+S=\epsilon\}.$ Making the following change of coordinates

we obtain the system

such that, under the coordinate change $k(y)=k(y+\varepsilon)$ with

Now, the splitting manifold becomes $\{y=0\}$. To make the computation easier, let us do a change of coordinates $(x,y)\rightarrow(y,x)$, which now the discontinuity manifold is given by $y-axis$. We will denote by $Z^{L}$ the vector field in the left side of the $y-axis$ and by $Z^{R}$ the vector field in the right side of the $y-axis$, where in our case, we obtain the following linear system:

The classical theory of linear differential systems shows that such types of systems do not have limit cycles. In this way, the only alternative for the existence of isolated periodic orbits is restricted to orbits that have some point in the discontinuity set.

As already observed, in planar piecewise smooth differential systems there are two kinds of limit cycles: sliding limit cycles and crossing limit cycles. Sliding limit cycles contain some segment of the discontinuity set, and crossing limit cycles only contain some points of the discontinuity set, which the associated orbit crosses the discontinuity set, in these points.

Here, we are interested only in crossing limit cycles. Then, in order to firsts analyzes, let us assume $\mathrm{x}=(x,y)^{T}\in\mathbb{R}^{2}$, $A^{L}=(a_{ij}^{L})$ and $A^{R}=(a_{ij}^{R})$ are $2\times 2$ constants matrices, and $b^{R}=(b_{1}^{R},b_{2}^{R})$ and $b^{L}=(b_{1}^{R},b_{2}^{L})$ are constant vectors of $\mathbb{R}^{2}$. Omitting the indexes $L,$ and $R$, we can write the linear piecewise system as

When we consider the discontinuity set being y-axis then we are considering the vector field with index $L$ for $x<0$ and the vector field with index $R$ for $x>0$. Analyzing the Lie derivative of each field on the y-axis, we see that the regions are determined by the product of $a_{12}^{L}y+b_{1}^{L}$ and $a_{12}^{R}y+b_{1}^{R}$.

Since our object of study is the (9) system we have $a_{12}^{L}=a_{12}^{R}=\alpha\,(1-\beta)$. If $\alpha\,(1-\beta)=0$, that is, $\alpha=0$ or $\beta=1$, then the first coordinate, of the vectors fields, has always the same sign on the line x=0, making it impossible to create limit cycles crossing the discontinuity line. Furthermore, simple calculations show us that the linear systems in (9) can only have an equilibrium of attractor node type, remember that $\beta>0>-k_{1}$. In summary, we have just proved the following proposition.

Assuming $a_{12}^{L}=a_{12}^{R}=\alpha\,(1-\beta)\neq 0$, we can regard another change of coordinates given by [12, 13], where is possible to write the linear system in terms of the determinant and the trace of the matrix, and which is described in the next proposition.

Omitting the indexes $L,$ and $R$, the systems in Proposition 4 can be written as

Now consider $[tr(A)]^{2}>det(A)>0$ and denote $\Delta=[tr(A)]^{2}-det(A)>0$. Then, considering $\lambda_{1}$, $\lambda_{2}$ the eigenvalues of the matrix $A$ with $\xi_{1}$, $\xi_{2}$ the associated eigenvectors then we can have

such that, clearly, $\lambda_{1}>\lambda_{2}.$ The equilibrium point and its invariant manifolds are

In this way, for $\alpha\,(1-\beta)\neq 0$, using the Proposition 4 we can write (9) as the following

such that $\Sigma=\{x=0\}$ and

Clearly, the vector field on the left $Z^{L}$ has the eigenvalues $\lambda_{1}^{L}=-k_{0}-\beta<0$ and $\lambda_{2}^{L}=-k_{0}-1<0$. On the other hand, the vector field on the right $Z^{R}$ has the eigenvalues $\lambda_{1}^{R}=-k_{1}-\beta<0$ and $\lambda_{2}^{R}=-k_{1}-1<0$. In the sequence, we will consider

where $\alpha^{R}_{\beta,\varepsilon}>\alpha^{L}_{\beta,\varepsilon}$, $\alpha^{R}_{\beta,\varepsilon}=\alpha^{L}_{\beta,\varepsilon}$ and $\alpha^{R}_{\beta,\varepsilon}<\alpha^{L}_{\beta,\varepsilon}$, respectively, if $\varepsilon<\varepsilon^{*}_{\beta}$, $\varepsilon=\varepsilon^{*}_{\beta}$ and $\varepsilon>\varepsilon^{*}_{\beta}$. Then, analyzing by direct calculations the equilibrium points, the Lie derivative of the system linear (12), and some analysis of signs, we have the following results.

The Proposition 6 assures us that it is not possible for the appearance of a fold of order greater than two for the system (12). This is because for $\alpha=\alpha^{L}_{\beta,\varepsilon}$ or $\alpha=\alpha^{R}_{\beta,\varepsilon}$, values which nullify $(Z^{\pm})^{2}f$ at the respective bending point, cause a collapse between of equilibrium point and the fold point in a boundary equilibrium of the respective system. And so, $(Z^{\pm})^{n}f=0$ for all $n$, since we have an equilibrium of the system $Z^{\pm}$ on the discontinuity line, that is a boundary equilibrium. Moreover, the invariant manifolds (15) and (17) intersect the discontinuity set $\Sigma=\{x=0\}$, respectively, at the points $(0,y_{m_{1}}^{L,R})$ and $(0,y_{m_{2}}^{L,R}+(k_{0}-k_{1})\,\varepsilon)$, where

Now we will introduce a new auxiliary function that will help us to study the periodic orbits of the piecewise system given by Proposition 4:

where $t$ is variable and $r_{1},r_{2}$ are parameters. Then we have the following lemma with proof elementary.

Turning our attention to the linear system (10) then its solution starting at $\mathrm{x}=(x_{0},y_{0})^{T}$, is given by

where $\lambda_{1}$ and $\lambda_{2}$ correspond to eigenvalues to the matrix $A$ such that $\lambda_{1}>\lambda_{2}$ and $(tr(A))^{2}>det(A)>0.$ And, in order to apply the Lemma 8 in the sequence, we need to fix a order relationship between the eigenvalues of Welander’s system (12).

• (I)

  The case $0<\beta<1$ guarantees that $0>\lambda_{1}^{L}>\lambda_{2}^{L}$ and $0>\lambda_{1}^{R}>\lambda_{2}^{R}$. Then we consider the auxiliary function (19), in the respective region, as

  (21)			$\displaystyle\psi_{\lambda_{1}^{L},\lambda_{2}^{L}}^{L}(t)=1-\beta-(k_{0}+1)\,e^{-(k_{0}+\beta)t}+(k_{0}+\beta)\,e^{-(k_{0}+1)t},\;\;\mbox{for}\;Z^{L},$
  			$\displaystyle\psi_{\lambda_{1}^{R},\lambda_{2}^{R}}^{R}(t)=1-\beta-(k_{1}+1)\,e^{-(k_{1}+\beta)t}+(k_{1}+\beta)\,e^{-(k_{1}+1)t},\;\;\mbox{for}\;Z^{R}.$

• (II)

  In turn, when $1<\beta$ we have $\lambda_{1}^{L}<\lambda_{2}^{L}<0$ and $\lambda_{1}^{R}<\lambda_{2}^{R}<0$. Then we consider the auxiliary function (19), in each region, as

  (22)			$\displaystyle\psi_{\lambda_{2}^{L},\lambda_{1}^{L}}^{L}(t)=\beta-1-(k_{0}+\beta)\,e^{-(k_{0}+1)t}+(k_{0}+1)\,e^{-(k_{0}+\beta)t},\;\;\mbox{for}\;Z^{L},$
  			$\displaystyle\psi_{\lambda_{2}^{R},\lambda_{1}^{R}}^{R}(t)=\beta-1-(k_{1}+\beta)\,e^{-(k_{1}+1)t}+(k_{1}+1)\,e^{-(k_{1}+\beta)t},\;\;\mbox{for}\;Z^{R}.$