Discontinuous differential equation for modelling the Antarctic Circumpolar Current

1.   Introduction

In recent years, the mathematical study of geophysical flows has been of great interest and has attracted much attention, particularly the equatorial flow, Antarctic Circumpolar Current (ACC), and Arctic gyres. We recommend readers to refer to monographs and the latest mathematical achievements to understand these oceanic phenomena. In addition, fluid stratification is an inherent characteristic of physical oceanography and is closely related to large-scale oceanic movements. The baroclinic instability is the main cause of the lateral and vertical stratification in ACC ^[1,2]. Moreover, the generation of stratification phenomenon also relies on eddies. This reason is that the interaction of imposed patterns of mechanical and buoyancy forcing by the eddies can establish stratification in both the horizontal and vertical directions ^[1,3].

For the motion of ACC, the existence of stratification that accommodates the observed sharp changes in water density (due to variations in temperature and salinity ^[4,5,6,7]), Constantin and Johnson ^[8,9] have constructed exact solutions to the governing equations of geophysical fluid dynamics in terms of spherical coordinates. These solutions represent purely azimuthal and depth-varying flows, which can be used to model equatorial undercurrent and ACC, respectively, and describe a purely homogeneous flow without stratification. For the investigation of stratified flow, Henry and Martin ^[10,11,12] constructed the exact solutions for geophysical fluid dynamics, which allows for the continuous stratification of equatorial flow that varies with the change of depth and latitude. Further, this method is extended to discontinuous density ^[13] and discontinuous varying density together with forcing terms in ACC ^[14,15].

In this paper, we investigate the existence of solutions for the stratified fluid flow corresponding to the different nonlinear vorticity in ACC. In response to the effects of density variation and forcing terms, we propose a new approach to determine the existence of stratified flow solutions. This method utilizes the Green's function and different nonlinear vorticity terms to derive the specific expression of the solution on both sides of the stratification and further proves the existence of stratification caused by density in the ACC.

The arrangement of this paper is presented as follows. In Section 2, we will state the governing equations with boundary conditions for ACC in terms of spherical coordinates. In Sections 3 and 4, based on the boundary conditions, we establish the existence of solutions corresponding to different nonlinear vorticity terms including piecewise constant coefficient, piecewise linear, and piecewise nonlinear, and verify the results graphically to represent the stratified flow in ACC. Section 5 extends to a general discontinuous case.

2.   Preliminary

Consider the spherical coordinates first, with $\theta$ $\in$ $[0, \pi)$ as the polar angle (with $\theta = 0$ corresponding to the North Pole) and $\varphi$ $\in$ $[0, 2\pi)$ as the angle of longitude (or azimuthal angle). In terms of the stream function $\psi(\varphi, \theta)$, the horizontal gyre flow on the spherical Earth has azimuthal and polar velocity components given by

Then, we introduce the stream function $\Psi (\theta, \varphi)$ associated with the vorticity of the ocean motion (the effects of the Earth's rotation are not taken into account), defined as

and the vorticity equation of the gyre flow is

where $\omega > 0$ is the nondimensional form of Coriolis parameter. $F(\Psi-\omega\cos\theta)$ and $2\omega \cos\theta$ are the oceanic vorticity and the planetary vorticity. By applying the stereographic projection of the unit sphere centered at the origin from the North Pole to the equatorial plane, the model (2.1) in spherical coordinates can be transformed into an equivalent plane elliptic partial differential equation, defined by

where $(r, \phi)$ are the polar coordinates in the equatorial plane. Thus (2.2) can be transformed into

Using the Cartesian coordinates $(x, y)$ in the complex $\xi$-plane, the above equation is equivalent the semi-linear elliptic partial differential equation

If the gyre flow has no variations in the azimuth, then the radially symmetric solutions of (2.3) have a form of $\psi = \psi(r)$. In terms of the change of variables $\psi = U(s)$, $s_1 < s < s_2$ with

for $0 < r_- < r_+ < 1$, the equation (2.3) can be transformed into the second-order differential equation

Considering the change of variables

we have

where $a(\cdot)$, $b(\cdot)$ are two positive continuous functions given by

and with the new boundary conditions

This means that the ACC is a streamline at the boundary. Moreover, we consider when $F(u)$ is a discontinuous function (including three cases), which means the ACC is stratified.

3.   The case $m < {\gamma} < n$

We consider in this section that

3.1.   Preliminary results

We first solve a simple discontinuous problem

where

for $m, n, {\gamma}\in{\mathbb R}$, and $h_\pm\in C[0, 1]$. We look for a $C^1-$smooth solution of (3.2) such that

${\bf{p1)}}$ There is an $t_1\in(0, 1)$ with $u(t_1) = {\gamma}$.

${\bf{p2)}}$ $u(t) < {\gamma}$ for $t\in(0, t_1)$ and ${\gamma} < u(t)$ for $t\in(t_1, 1)$.

We see that if such a $u(t)$ exists, then it must have by p1) a form

where $G_\pm(t_1, t, s)$ are Green functions

and

Clearly,

A condition $u'_-(t_{1}) = u'_+(t_{1})$, $t_1\in(0, 1)$ is equivalent to

Hence, we solve

Since

there is a $t_1\in(0, 1)$ solving (3.5). We take such $t_1$. To satisfy p2), we suppose either

or

If (3.6) holds, then $u_+(t)$ is strictly concave on $[t_1, 1]$, which implies $u_+(t) > n+\frac{{\gamma}-n}{1-t_1}(1-t) > {\gamma}$ for $t\in(t_1, 1]$, and $u_+(t)$ has at most one extreme point in $[t_1, 1]$. Hence, $u_+'(t_1) > 0$ and, thus, $u_-'(t_1) > 0$. Hence, $u_-(t) < {\gamma}$ for $t < t_1$ near $t_1$. If there is an $t_2\in(0, t_1)$ with $u_-(t_2) = {\gamma}$, then there are $t_3\in(0, t_2]$ and $t_4\in(t_2, t_1)$ with $u_-'(t_3) = u_-'(t_4) = 0$, so there is $t_5\in(t_3, t_4)$ with $u''(t_5) = 0$, which contradicts to (3.6). Thus, $u_-(t) < {\gamma}$ for $0\le t < t_1$. Similar arguments work under (3.7). Summarizing, we arrive at the following result.

Theorem 3.1. Under assumption (3.1), either (3.6) or (3.7) holds, then there is a solution (3.3) of (3.2) with properties of p1) and p2), where $t_1$ solves (3.5).

3.2.   Piecewise constant oceanic vorticity

We apply results from Subsection 3.1 by considering

where $b_1, b_2, \gamma$ are constants. Then,

where

Since $-\frac{c(t)}{a(t)} = 2 \frac{1-e^{t(s_2-s_1)+s_1}}{1+e^{t(s_2-s_1)+s_1}}$, set

Now, condition (3.6) means

and condition (3.7) means

By Theorem 3.1, we get the following result.

Theorem 3.2. Consider (3.8) and assumption (3.1). If either (3.10) or (3.11) holds, then there is a solution (3.3) of (3.2) with properties of p1) and p2), where $t_1$ solves (3.5).

Note that now (3.5) is a linear function of $b_1$, $b_2$, ${\omega}$, $m$, $n$, and ${\gamma}$ when $t_1\in(0, 1)$ is considered as a parameter of the form

So, we get a surface in the space $(t_1, b_1, b_2, {\omega}, m, n)$ when either (3.10) or (3.11) holds. The righthand side of (3.12) has an explicit formula, but it is rather awkward, so we do not go into more details, we just add a concrete example.

Example 3.3. We take

Clearly (3.11) holds for these values (3.13), since $p_1 = 2\frac{1-e^{20}}{1+e^{20}} < 0$ and $p_2 = \frac{1-e^{1}}{1+e^{1}} < 0$. We solve ${\gamma} = {\gamma}(t_1)$ from (3.12) with (3.13) (see Figure 1), which is again a huge formula not suitable for writing it here. However, we numerically find that ${\gamma}(t_1)$ is increasing on $[0, 0.0275]$ from $1$ to $2.72678$, then decreasing on $[0.0275, 0.101948]$ from $2.72678$ to $0.538563$, and increasing on $[0.101948, 1]$ from $0.538563$ to $150$. This gives that for ${\gamma}\in(0, 2.72678)$, we have $3$ solutions of (3.12), then $2$ solutions for ${\gamma} = 2.72678$ and $1$ solution for ${\gamma}\in(2.72678,150)$. Thus, ${\gamma}$ has a cusp type bifurcation.

3.3.   Piecewise linear oceanic vorticity

Now, we extend results of Subsection 3.2 by considering a piecewise linear vorticity case

where $k_1$, $k_2$, $b_1$, $b_2$ are constants. Then,

We take $\mu < m < n < \nu$ and consider a subset

Set

Then, we derive

for any $u\in{\Omega}_1$. Now, condition (3.6) means

and condition (3.7) means

For any $t_1\in[0, 1]$ and $u\in {\Omega}_1$, we consider (3.3) with (3.15). Then, (3.4) and (3.18) imply

By setting

(3.21) gives

By assuming

(3.23) implies

We denote (3.3) with (3.15) by $\Psi(u, t_1)(t)$, i.e,

By following (3.5), we set $\psi : {\Omega}\to{\mathbb R}$, ${\Omega} = {\Omega}_1\times(0, 1)$ as

see (3.15). Summarizing, we obtain the following result.

Theorem 3.4. Assuming either (3.19) or (3.20), and (3.24) with (3.22), there is a solution (3.3) of (3.2) with properties of p1), p2), and (3.16).

Proof. It is standard to verify that

is continuous and compact. We introduce a homotopy $H = (H_1, H_2): [0, 1]\times{\Omega}\to C_0\times{\mathbb R}$ as follows:

Clearly, (3.25) gives

and noting

we see that

for a small ${\delta} > 0$ and ${\Omega}_{\delta} = {\Omega}_1\times ({\delta}, 1-{\delta})$. This means that we can compute a Leray-Schauder degree

Thus, there is an $(u, t_1)\in{\Omega}_{\delta}$ such that

By (3.19), (3.20), and Theorem 3.1, we get a desired solution. The proof is finished.

Next, it is clear that $h_+$ and $h_-$ in (3.15) are Lipschitz continuous in $u$ with constants $|k_1|$ and $|k_2|$, respectively. By introducing

we can easily check that $\Psi(u, t_1)$ is Lipschitz continuous in $u$ on ${\Omega}_1$ with a constant

uniformly with respect to $t_1\in(0, 1)$. By using a Banach fixed point theorem, we can improve Theorem 3.4 as follows.

Theorem 3.5. Assume either (3.19) or (3.20), and (3.24) with (3.22). If

then for any $t_1\in(0, 1)$, there is a unique fixed point of

with (3.16). Inserting it into (3.27), we get a scalar equation,

It has a solution which can be approximately located by

Thus, a solution (3.3) of (3.2) with properties of p1) and p2) can be approximately located.

Example 3.6. We extend Example 3.3 for (3.14), satisfying additional conditions. To begin, we suppose (3.30) with (3.29). Then, (3.26) is a global contraction on $C_0$ and we split it to

for

and

We solve (3.31) by the Banach fixed point theorem to get

Finally, we get (3.32), which is solvable on $(0, 1)$ according to the above arguments. Clearly, (3.36) gives

Next, (3.13) implies

thus applying also (3.15) and (3.37), we obtain

if

Summarizing, we arrive at the following conclusion.

Theorem 3.7. Consider (3.14) with parameters from Example 3.3. Suppose (3.30) and (3.38) with (3.29). Then, (3.2) has a solution with properties of p1) and p2).

3.4.   Piecewise nonlinear oceanic vorticity

Finally, we extend results of Subsection 3.3 by considering a nonlinear case given as

where $F_1$, $F_2$ are continuous. Then,

We again consider a subset $\Omega_1$ from (3.16). We set

Clearly, (3.26) now takes a form

Now we can directly follow step by step arguments of Subsection 3.3 to extend the above results for a nonlinear case. So, we get the following results.

Theorem 3.8. Theorem 3.4 is valid for (3.39) with (3.40).

Theorem 3.9. Assuming

Theorem 3.5 is valid for (3.39) with (3.40).

Example 3.10. We extend Example 3.6 for

To begin, we suppose (3.30) with (3.42). Then, (3.41) is a global contraction on $C_0$ and we can directly repeat the argument of Example 3.6 to get the following result.

Theorem 3.11. Consider (3.39) with parameters from Example 3.3. Suppose (3.30) with (3.42) and also (see (3.38))

Then, (3.2) has a solution with properties of p1) and p2).

By using (3.13), we compute from (3.28)

and (3.29) has a value

so (3.30) reads

Furthermore, we calculate

Hence, (3.44) becomes

The smallness of the above values is caused by the exponential terms in $a(t)$, $c(t)$, and $s_2 = 20$.

Remark 3.12. A reverse inequality of (3.1) given by

is symmetric to (3.1) by taking a symmetric transformation $t\leftrightarrow 1-t$, $t\in[0, 1]$. So considering $u(t)\leftrightarrow u(1-t)$, results for the case (3.1) in this section are extended straightforwardly to the case (3.45).

4.   The case $m < {\gamma}$ and $n < {\gamma}$

We consider in this section that

This situation (4.1) is different from both (3.1) and (3.45), but we can extend a method of Section 3.

4.1.   Preliminary results

We look for a $C^1-$smooth solution of (3.2) such that

${\bf{p3)}}$ There are $0 < t_1 < t_2 < 1$ with $u(t_1) = u(t_2) = {\gamma}$.

${\bf{p4)}}$ $u(t) < {\gamma}$ for $t\in(0, t_1)\bigcup(t_2, 1)$ and ${\gamma} < u(t)$ for $t\in(t_1, t_2)$.

We see that if such a $u(t)$ exists, then it must have by p3) a form

where $G_\pm$ are Green functions

and

Clearly, $G_\pm < 0$. A condition $u'_-(t_{1}) = u'_+(t_{1})$, $0 < t_1 < t_2 < 1$ is equivalent to

Next, a condition $u'_-(t_{2}) = u'_+(t_{2})$, $0 < t_1 < t_2 < 1$ is equivalent to

We consider

Since

we derive

By supposing (3.6) and following arguments in Subsection 3.1, we can see that p4) is satisfied.

Of course, solving (4.3) and (4.4) is not easy. If there is an $t_1\in(0, 1)$ such that

then using (4.3), (4.4), (4.5), and (3.6), we see that $h_-(t) < 0$ on $[0, 1]$ and $u_-'(t_1-) = u_-'(t_1+)$. This gives a solution (4.2) for $t_1 = t_2$, which is grazing the line $u = {\gamma}$ at $t = t_1$.

On the other hand, if

which holds if

then we can consider a Brouwer degree or a winding number of a planar vector field $(\psi_1, \psi_2)$

which can be computed by using (4.5). Summarizing, we have the following.

Theorem 4.1. Under assumption (4.1), (3.6) holds, then there is a solution (3.3) of (3.2) with properties of p3) and p4), where $t_1$ and $t_2$ solves (4.3) and (4.4). Moreover, such $t_1$ and $t_2$ exist if supposing (4.7) and

holds.

Next, we compute (see (4.5))

along a counterclockwise oriented border

We see that the vector field $\{(\psi_1(0, t), \psi_2(0, t))\}_{t\in[0, 1]}$ is located in the right half plane $x_1 > 0$, while the vector field $\{(\psi_1(t, 1), \psi_2(t, 1))\}_{t\in[0, 1]}$ is located in the lower half plane $x_2 < 0$. Thus, there is no whole rotation from a vector

to a vector

along the borders $\{[0, t], [t, 1]\mid t\in[0, 1]\}$. As a matter of fact, we can use along the borders $\{[0, t], [t, 1]\mid t\in[0, 1]\}$ in the following homotopy:

which is suitable, since $h_{l1}(t, {\lambda}) = ({\gamma}-m)\ne0$ and $h_{r2}(t, {\lambda}) = (n-{\gamma})\ne0$. Note that $h_{l}(t, 1)$, $h_{r}(t, 1)$ gives linear lines connecting $\psi(0, 0)$ with $\psi(1, 1)$ via $({\gamma}-m, n-{\gamma})$.

Next, under (4.8), the vector field $\{(\psi_1(t, t), \psi_2(t, t))\}_{t\in[0, 1]}$ is located into a quadrant $x_1 > 0$, $x_2 < 0$, which gives that there is no whole rotation along a counterclockwise oriented $\partial\triangle$. This means

Consequently, under assumptions

we cannot use a topological degree argument of Theorem 4.1, but Theorem 4.1 is still applicable. We introduce a mapping

and consider the set

Note $\chi(t_1, t_2) = \psi(t_1, t_2)-({\gamma}-m, n-{\gamma})$. So for $(t_1, t_2)\in\Upsilon$, we take

To complete our discussion, we assume

Then, (4.10) shows that both $\psi_{1, 2}(t, t)$ are decreasing on $[0, t]$. So, if (4.7) holds along with

then the vector field $\psi$ along a counterclockwise oriented border $\partial\triangle$ takes one counterclockwise turn, so (4.9) holds with $= 1$. Otherwise, (4.11) holds. Furthermore, (4.16) requires $\psi(1, 1) < 0$, i.e.,

Since $\psi_1(0, 0) = {\gamma}-m > 0$, under (4.17), there is a unique $t_0\in(0, 1)$ such that $\psi_1(t_0, t_0) = 0$, i.e.,

Thus, $\psi_1(t, t) < 0$ for $t\in(t_0, 1]$ and $\psi_1(t, t) > 0$ for $t\in[0, t_0)$. If $\psi_2(t_0, t_0)\le0$, then $\psi_2(t, t)\le0$ for $t\in(t_0, 1]$. Hence, (4.16) is not satisfied. On the other hand, if $0 < \psi_2(t_0, t_0)$, i.e.,

then (4.16) holds. Clearly, (4.7) also holds. Summarizing, we arrive at the following result.

Lemma 4.2. Under assumptions (4.15), (4.17), and (4.19), where $t_0$ uniquely solves (4.18), we have

4.2.   Piecewise constant oceanic vorticity

We apply results from Subsection 4.1 by considering (3.8). We can follow the method of Subsection 3.2.

Theorem 4.3. Consider (3.8) and assumption (4.1). If (3.10) holds, there is a solution (3.3) of (3.2) with properties of p3) and p4), when $t_1$ and $t_2$ solves (4.3) and (4.4).

Now, (4.3) and (4.4) have forms

and

Again, (4.20) and (4.21) are linear functions of $b_1$, $b_2$, ${\omega}$, $m$, $n$, and ${\gamma}$ when $(t_1, t_2)\in\triangle$ are considered as parameters. So, we get surfaces in the space $(t_1, t_2, b_1, b_2, {\omega}, m, n)$ when (3.10) holds. The righthand sides of (4.20) and (4.21) have explicit formulae, but they are rather awkward, so we do not go into more details, but just add a concrete example.

Example 4.4. We take

Then

and (3.10) gives (4.15). We recall (4.14). We numerically check that $\chi_1(t, t)$ is decreasing on $[0, 1]$ from $0$ to $-0.579183$, while $\chi_2(t, t)$ is decreasing on $[0, 1]$ from $26.643$ to $0$ (see Figures 2 and 3). This means $(t, t)\in\Upsilon$, $t\in[0, 1]$ and that for a suitable $({\gamma}-m, n-{\gamma})$, we have (4.9) with $= 1$: $-({\gamma}-m, n-{\gamma})$ from a region between the graph of $\chi(t, t)$ and left axis $x_1$ and upper axis $x_2$, so (4.16) holds, for instance, we can take $({\gamma}-m, n-{\gamma}) = (0.2, -0.5)$. On the other hand, we see from Figure 3 that there are $({\gamma}-m, n-{\gamma})$ when (4.11) holds, but Theorem 4.1 is still applicable, i.e, (4.14) holds for some $(t_1, t_2)\in\Upsilon$: $-({\gamma}-m, n-{\gamma})$ from a region $\chi(\Upsilon)$ on the left side of the graph of $\chi(t, t)$, for instance, we can take $({\gamma}-m, n-{\gamma}) = (0.7, -0.5)$. These results are consistent with our discussion at the end of Subsection 4.1.

4.3.   Piecewise nonlinear oceanic vorticity

Finally, we extend results of Subsection 3.4 by considering a nonlinear case given as (3.39). From (4.2), (3.41), now possesses a form

We take

We again consider a subset $\Omega_1$ from (3.16) and (3.40). By setting

similar estimates like (3.21) give

By assuming

(4.25) implies

Note, (4.23) can be continuously extended to $0 < t_1 = t_2 < 1$. Summarizing, we obtain the following result.

Theorem 4.5. Assuming (4.26) with (4.24), there is a fixed point $u(t_1, t_2) = \Psi(u(t_1, t_2), t_2, t_1)$ in (3.16) for any $(t_1, t_2)\in\triangle_D = \{0 < t_1\le t_2 < 1\}$.

Proof. It is standard to verify that

is continuous and compact. Clearly, (4.27) gives

The proof follows from the Schauder fixed point theorem.

Next, by introducing

and assuming (3.42), we can easily check that $\Psi(u, t_2, t_1)$ is Lipschitz continuous in $u$ on ${\Omega}_1$ with a constant (3.29) uniformly with respect to $t_1, t_2\in\triangle_D$. By using a Banach fixed point theorem, Theorem 3.5 states as follows.

Theorem 4.6. If (3.30) holds from (3.29), with (4.28) in addition in Theorem 4.5, then $u(t_1, t_2)$ is unique.

Inserting $u(t_1, t_2)$ into (4.3) and (4.4), we get equations

and

We arrive at the following result.

Theorem 4.7. Assume (3.19) in addition in Theorems 4.5 and 4.6. If $0 < t_1 < t_2 < 1$ solves (4.29) and (4.30), then $u(t_1, t_2)$ solves (3.2) with (4.1) for (3.15), and it has properties of p3) and p4).

Solving (4.29) and (4.30) is difficult. We can use an approximation method of Theorem 3.5 for Theorem 4.6 as follows:

Furthermore, following arguments of Subsection 4.1, we get

for $t\in(0, 1)$. To get (4.15), we assume (see (3.19))

At the end, we have to compute

for large $n$ by using a method from Subsection 4.1.

Remark 4.8. We can take a concrete example by extending Example 4.4 for $F_1(u) = k_1|\sin u| +b_1$ and $F_2(u) = k_2|\cos u|+b_2$ with sufficiently small $k_1 > 0$ and $k_2 > 0$. However, we do not present details since they are similar.

Remark 4.9. A reverse inequality of (4.1) given by

leads to (4.1) by taking a transformation $u\leftrightarrow -u$, and the results of this section can be directly reformulated.

5.   Final remarks

In the future work, we can consider general discontinuities, then we must pass to differential inclusions via Filippov convexification of (3.39)

where

We rewrite (5.1) as a fixed point problem

where

and

So, we introduce a multifunction

and (5.2) is equivalent to

Now, if $F_1(u)$ and $F_2(u)$ are globally bounded and continuous, then the Kakutani fixed point theorem gives the existence of the solution of (5.3) ^[16,17,18]. The advantage of the above Sections 3 and 4 is finding detailed solutions, which can be constructed or approximated, while (5.3) gives a vague form of solutions. For a numerical analysis of (5.1), we use an approximation of (5.2) given by

for $\varepsilon > 0$ sufficiently small.

Author contributions

Michal Fečkan: Writing-Original draft preparation, Writing-Reviewing and Editing; Shan Li: Writing-Reviewing and Editing; JinRong Wang: Conceptualization, Writing-Reviewing and Editing.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work is partially supported by the National Natural Science Foundation of China (12371163), by the Slovak Research and Development Agency under the Contract no. APVV-23-0039, and by the Slovak Grant Agency VEGA No. 1/0084/23 and No. 2/0062/24.

Conflict of interest

The authors declare there is no conflict of interest.