# Opinion Dynamics Optimization Through Noncooperative Differential Games
††thanks: This work was supported by SGS, VŠB - Technical University of
Ostrava, Czech Republic, under grant No. SP2023/12 “Parallel processing of Big
Data X”.

Hossein B. Jond  Department of Computer Science
VŠB-Technical University of Ostrava
Ostrava-Poruba, Czech Republic
<EMAIL_ADDRESS>
###### Abstract

In this paper, I study optimizing the opinion formation of a social network of
a population of individuals on a graph whose opinion evolves according to the
Hegselmann-Krause model for opinion dynamics. I propose an optimization
problem based on a differential game for a population of individuals who are
not stubborn. The objective of each individual is to seek an optimal control
policy for her own opinion evolution by optimizing a personal performance
index. The Nash equilibrium actions and the associated opinion trajectory with
the equilibrium actions are derived for the opinion optimization model using
Pontryagin’s principle. The game strategies were executed on the well-known
Zachary’s Karate Club social network. The resulting opinion trajectories
associated with the game strategies showed that in non-stubborn Zachary’s
network, the opinions moved toward the average opinion of the network, but a
consensus of final opinions did not necessarily emerge.

###### Index Terms:

game theory, Hegselmann-Krause model, opinion dynamics, optimization

##  I Introduction

Opinion dynamics, with its origins in sociology, is the study of the dynamical
processes of public opinion formation, diffusion, and evolution [1, 2].
Opinion dynamics have rapidly spread beyond sociology into other disciplines
such as physics, mathematics, computer science, and control theory, to name a
few. Among the various proposed models, the bounded confidence models of
opinion dynamics have attracted more attention [3, 4]. The bounded confidence
concept was incorporated into opinion dynamics by Hegselmann and Krause [4].
In the Hegselmann-Krause (HK) model, a finite number of agents fuse their
opinions only with those of others whose opinions do not differ more than
their confidence bound. Several works have studied the extensions of the HK
model under various settings [7, 5, 6].

Opinion dynamics optimization is an emerging topic of interest for those
dealing with social networks. In [8], an optimization procedure is presented
to select informed agents to prevent the influence of the annoying agents in
the network. In another study, the problem of optimizing the placement of
stubborn agents in a social network with the aim of maximally influencing the
population was studied [9].

One interesting approach to studying opinion dynamics in social networks is a
game-theoretic approach. In [7], the authors showed the HK model can be
formulated as a sequence of the best response dynamics of a potential game.
The work [10] investigated a non-cooperative differential game model of
opinion dynamics with an open-loop information structure. The opinion dynamics
of a multiple-population social network were investigated through the
application of a multiple-population mean field game in [11]. The work [12]
investigated the spread of opinions in a binary opinion model as an
evolutionary game. A differential game in [13] was used to model the opinion
behavior of stubborn agents in a social network in the presence of a troll.

In this study, I model and evaluate the opinion behavior of selfish and self-
interested agents in a social network whose opinion fusion rules conform to
the HK opinion dynamics model and are constrained by a communication graph
using non-cooperative differential games. The simulation results indicate that
a network whose population has adopted game equilibrium strategies does not
have a consensus of final opinions.

The main contributions of this paper in comparison with similar works in [10,
13] are twofold. First, [10, 13] considered the single integrator dynamics
model for opinion evolution, whereas I utilize the HK opinion dynamics model.
The HK model, in contrast to the model with a single integrator, takes into
account the topological structure of the social network. Second, the cost
functions in [10, 13] involve minimizing disagreements in the network
throughout the entire opinion formation process. In my optimization, I take
final opinion disagreements into account. This choice of cost function results
in more concise explicit expressions for the Nash equilibrium and the opinion
trajectories associated with it. In addition, the distributed information
structure of Nash equilibrium is observed.

The paper is organized as follows. In Section II, I present a differential
game model of the HK opinion dynamics for opinion optimization. In Section
III, I derive the open-loop Nash equilibrium solution and the associated
opinion trajectory with it for the optimization model. In addition, I also
present an optimal control model of the HK opinion dynamics. In Section IV,
the results from the previous section are implemented on a real-world social
network to observe the evolution of individuals’ opinions. Conclusions and
future works are discussed in Section V.

##  II Opinion Dynamics Optimization

Consider a social network of $n$ agents indexed $1$ through $n$ on a
communication graph $\mathcal{G}(\mathcal{V},\mathcal{E})$. The set of
vertices $\mathcal{V}=\\{1,\cdots,n\\}$ corresponds to the set of agents. Each
edge $(i,j)\in\mathcal{E}$ represents a mutual opinion flow between node $i$
and node $j$. The set of neighbors of vertex $i$ is denoted by
$\mathcal{N}_{i}=\\{j\in\mathcal{V}:(i,j)\mbox{ or }(j,i)\in\mathcal{E},j\neq
i\\}$. I make the following assumption. The social graph $\mathcal{G}$ is
connected. The connectivity of $\mathcal{G}$ means at least one globally
reachable node (a root node of a spanning tree on the graph). The connected
social graph means that each agent has at least one neighbor with whom they
mutually interact (i.e., $\mathcal{N}_{i}\neq\emptyset$ or
$|\mathcal{N}_{i}|\neq 0$ for all $i\in\mathcal{V}$).

Let $x_{i}(t)\in\mathit{\mathbb{R}}$ be the opinion of agent $i$ at time
$t\in[0,t_{f}]$ where $t_{f}$ is a terminal time. In the HK model, the
evolution of $x_{i}(t)$ at each stage $k=0,1,\cdots$ is as follows

$x_{i}(k+1)=\frac{1}{|\mathcal{N}_{i}|}\sum_{j\in\mathcal{N}_{i}}x_{j}(k).$
(1)

In this model, each agent’s opinion at each stage is the average opinion of
her graph neighbors. The graph neighbors here are equivalent to the bounded
confidence concept in the original HK model. An agent with the confidence
bound $\epsilon_{i}$ takes only those agents $j$ into account whose opinions
differ from her own not bigger than her confidence bound. The set of such
agents $j$ is denoted by
$N_{i}(t)=\\{j\in\mathcal{V}:|x_{i}(t)-x_{j}(t)|\leq\epsilon_{i}\\}$ which in
terms of a time-invariant social graph $\mathcal{G}(\mathcal{V},\mathcal{E})$,
is equivalent to the set of the graph neighbors for agent $i$ as defined
before, i.e., $\mathcal{N}_{i}$. Another interpretation for $\mathcal{N}_{i}$
could be the set of agents whom a particular agent $i$ trusts to share and
fuse opinion [14].

Define $u_{i}(t)\in\mathit{\mathbb{R}}$ as agent $i$’s influence effort or
simply her control input at time $t\in[0,t_{f}]$. The HK model (1) with the
control input in continuous time is given by

$\displaystyle \dot{x}_{i}(t)$
$\displaystyle=\frac{1}{|\mathcal{N}_{i}|}\sum_{j\in\mathcal{N}_{i}}\big{(}x_{j}(t)-x_{i}(t)\big{)}+b_{i}u_{i}(t)$
$\displaystyle=\frac{1}{|\mathcal{N}_{i}|}\big{(}\sum_{j\in\mathcal{N}_{i}}x_{j}(t)-|\mathcal{N}_{i}|x_{i}(t)\big{)}+b_{i}u_{i}(t),$

or simply,

$
\dot{x}_{i}(t)=\frac{1}{|\mathcal{N}_{i}|}\sum_{j\in\mathcal{N}_{i}}x_{j}(t)-x_{i}(t)+b_{i}u_{i}(t),$
(2)

where $b_{i}$ is a nonzero constant.

The HK model assumes that the social network population is not stubborn, which
means that the individuals do not hold any prejudices. In other words, their
opinion evolves without any influence from their initial opinions, which suits
modeling a social network population whose attitudes change independently of
their previous beliefs. The process of opinion formation in a rational, non-
stubborn social network at the individual level can be subject to each agent
attempting to minimize her disagreement with her graph neighbors while, in the
meantime, expending the least amount of influence effort. An appropriate cost
function that characterizes such behavior or preferences in a social network
for agent $i\in\mathcal{V}$ to minimize is

$J_{i}=\frac{1}{|\mathcal{N}_{i}|}\sum_{j\in\mathcal{N}_{i}}\big{(}x_{i}(t_{f})-x_{j}(t_{f})\big{)}^{2}+\int_{0}^{t_{f}}r_{i}u_{i}^{2}(t)~{}\mathrm{dt},
$ (3)

where $r_{i}\in\mathit{\mathbb{R}}$ is a positive scalar ($r_{i}>0$). The cost
function in (3) has two terms. The first term averages the sum of
disagreements between the final opinions of each agent and her graph
neighbors. The second term is the weighted control or influence effort made
during the entire opinion formation process.

In the control community, the optimization problem that emerged in this work
is known as differential game problems [15]. In the context of a differential
game, each agent of the network is referred to as a player. In this context,
each player seeks the control $u_{i}(t)$ that minimizes her cost function
$J_{i}$ with the given initial opinions, subject to the continuous-time HK
opinion evolution equation (2). In other words, the players in the game seek
to minimize their cost functions in order to find their control or influence
strategies $u_{i}(t)$ while their opinions evolve according to the
differential equation (2). The behavior of self-interested players in the game
of opinion formation is best reflected via noncooperative game theory. Under
the framework of noncooperative games, the players can not make binding
agreements, and therefore, the solution (i.e., the Nash equilibrium) has to be
self-enforcing, meaning that once it is agreed upon, nobody has the incentive
to deviate from [16]. In the next section, I derive the open-loop Nash
equilibrium solution for the optimization problem in (2) and (3).

##  III Main Result

Nash equilibrium is the main solution concept in noncooperative game
scenarios. A Nash equilibrium is a strategy combination of all players in the
game with the property that no one can gain a lower cost by unilaterally
deviating from it. The open-loop Nash equilibrium is defined as a set of
admissible actions ($u_{1}^{*},\cdots,u_{n}^{*}$) if for all admissible
($u_{1},\cdots,u_{n}$) the inequalities
$J_{i}(u_{1}^{*},\cdots,u_{i-1}^{*},u_{i}^{*},u_{i+1}^{*},\cdots,u_{n}^{*})\leq
J_{i}(u_{1}^{*},\cdots,u_{i-1}^{*},u_{i},u_{i+1}^{*},\cdots,u_{n}^{*})$ hold
for $i\in\\{1,\cdots,n\\}$ where $u_{i}\in\Gamma_{i}$ and $\Gamma_{i}$ is the
admissible strategy set for player $i$. The noncooperative differential game
and the unique Nash equilibrium associated with it are discussed in [15].

In the following, I present the main result, which is the open-loop Nash
equilibrium solutions and the associated opinion trajectories with the
equilibrium actions for the previously introduced optimization problem. Before
that, I define the following vectors and matrices to restate the optimization
problem in a compact form.

Define $A_{i}=[a_{ij}]$ where

$a_{ij}=\left\\{\begin{array}[]{ll}1&\mbox{if }(i,j)\in\mathcal{E}\mbox{ and
}i\neq j,\\\ 0&\mbox{if }(i,j)\notin\mathcal{E}\mbox{ and }i\neq j,\\\
0&\mbox{if }i=j,\end{array}\right.$

to be the adjacency matrix for each agent $i\in\mathcal{V}$. The degree matrix
is $D_{i}=\mathrm{diag}(0,\cdots,|\mathcal{N}_{i}|,\cdots,0)$ where
”$\mathrm{diag}{\textperiodcentered}$” stands for diagonal matrix. The graph
Laplacian matrix for each agent $i$ is defined as

$L_{i}=D_{i}-A_{i}.$ (4)

The global adjacency, degree, and Laplacian matrices are $A=\sum_{i}A_{i}$,
$D=\sum_{i}D_{i}$, and $L=\sum_{i}L_{i}$. All the aforementioned matrices are
symmetric. Additionally, let
$W_{i}=\mathrm{diag}(0,\cdots,\omega_{ii},\cdots,0)$ and $W=\sum_{i}W_{i}$.
Define vectors $x(t)=[x_{1}(t),\cdots,x_{n}(t)]^{\top}$,
$B_{i}=[0,\cdots,b_{i},\cdots,0]^{\top}$, and matrix $\Lambda=D^{-1}A-I$.

The optimization in (2) and (3) is restated in compact form as follows

$\displaystyle\min_{u_{i}}~{}J_{i}(u_{i})=\frac{1}{|\mathcal{N}_{i}|}x^{\top}(t_{f})L_{i}x(t_{f})+\int_{0}^{t_{f}}r_{i}u_{i}^{2}(t)~{}\mathrm{dt},$
(5) $\displaystyle{s.t.}$ $\displaystyle\dot{x}(t)=\Lambda
x(t)+\sum_{i=1}^{n}B_{i}u_{i}(t),\quad x(0)=x_{0},\quad i=1,\ldots,n.$

The presence of the Laplacian $L_{i}$ in the cost function above is due to its
sum-of-squares property (see [17]).

$
x^{\top}(t_{f})L_{i}x(t_{f})=\sum_{j\in\mathcal{N}_{i}}\big{(}x_{i}(t_{f})-x_{j}(t_{f})\big{)}^{2}.$

######  Theorem 1.

The unique Nash equilibrium actions and the associated opinion trajectory with
these equilibrium actions for opinion formation of a social network as the
noncooperative differential game in (5) are given by

$\displaystyle u_{i}(t$
$\displaystyle)=-\frac{1}{r_{i}|\mathcal{N}_{i}|}B_{i}^{\top}\mathrm{e}^{(t_{f}-t)\Lambda^{\top}}L_{i}H^{-1}(t_{f})\mathrm{e}^{t_{f}\Lambda}x_{0},$
(6) $\displaystyle x(t)$
$\displaystyle=\Big{(}\mathrm{e}^{t\Lambda}-\Psi(t)\Delta
H^{-1}(t_{f})\mathrm{e}^{t_{f}\Lambda}\Big{)}x_{0},$ (7)

where

$\displaystyle H(t_{f})=I+\Psi(t_{f})\Delta,$ (8)
$\displaystyle\Psi(t)=[\Psi_{1}(t),\cdots,\Psi_{n}(t)],$ (9)
$\displaystyle\Delta=[\frac{1}{|\mathcal{N}_{1}|}L_{1},\cdots,\frac{1}{|\mathcal{N}_{n}|}L_{n}]^{\top},$
(10)
$\displaystyle\Psi_{i}(t)=\int_{0}^{t}\mathrm{e}^{(t-\tau)\Lambda}S_{i}\mathrm{e}^{(t-\tau)\Lambda^{\top}}~{}\mathrm{d\tau},\quad
S_{i}=\frac{1}{r_{i}}B_{i}B_{i}^{\top}.$ (11)

###### Proof.

Define the Hamiltonian

$\mathcal{H}_{i}=r_{i}u_{i}^{2}(t)+\lambda_{i}^{\top}(t)\Big{(}\Lambda
x(t)+\sum_{i=1}^{n}B_{i}u_{i}(t)\Big{)},$ (12)

where $\lambda_{i}(t)$ is a co-state vector. According to Pontryagin’s
principle, the necessary conditions for optimality are
$\frac{\partial\mathcal{H}_{i}}{\partial u_{i}}=0$ and
$\dot{\lambda}_{i}(t)=-\frac{\partial\mathcal{H}_{i}}{\partial x}$. Applying
the necessary conditions on the Hamiltonian yield

$\displaystyle u_{i}(t)=-\frac{1}{r_{i}}B_{i}^{\top}\lambda_{i}(t),$ (13)
$\displaystyle\dot{\lambda}_{i}(t)=-\Lambda^{\top}\lambda_{i}(t),$ (14)

with the terminal condition

$\lambda_{i}(t_{f})=\frac{1}{|\mathcal{N}_{i}|}L_{i}x(t_{f}).$ (15)

The solution of (14) is uniquely determined by

$\lambda_{i}(t)=\mathrm{e}^{(t_{f}-t)\Lambda^{\top}}\lambda_{i}(t_{f}).$ (16)

Substituting this solution in (13) yields

$u_{i}(t)=-\frac{1}{r_{i}}B_{i}^{\top}\mathrm{e}^{(t_{f}-t)\Lambda^{\top}}\lambda_{i}(t_{f}).$
(17)

Substituting (13) into the opinion dynamics in (5) and then using (17), I get

$\displaystyle\dot{x}(t)$ $\displaystyle=\Lambda
x(t)-\sum_{i=1}^{n}S_{i}\lambda_{i}(t)$ $\displaystyle=\Lambda
x(t)-\sum_{i=1}^{n}S_{i}\mathrm{e}^{(t_{f}-t)\Lambda^{\top}}\lambda_{i}(t_{f}).$
(18)

The solution of (III) at $t$ is given by

$x(t)=\mathrm{e}^{t\Lambda}x_{0}-\sum_{i=1}^{n}\Psi_{i}(t)\lambda_{i}(t_{f}),$
(19)

where $\Psi_{i}(t)$ is defined in (11). Using notation (9), equation (19) is
rewritten as

$x(t)=\mathrm{e}^{t\Lambda}x_{0}-\Psi(t)\lambda(t_{f}).$ (20)

Stacking (15) for $i=1,\cdots,n$ yields

$\lambda(t_{f})=\Delta x(t_{f})$ (21)

with $\Delta$ defined in (10) and
$\lambda(t_{f})=[\lambda_{1}^{\top}(t_{f}),\cdots,\lambda_{n}^{\top}(t_{f})]^{\top}$.
Substituting $\lambda(t_{f})$ from (21) then into (20) at $t_{f}$ yields

$\displaystyle x(t_{f})=\mathrm{e}^{t_{f}\Lambda}x_{0}-\Psi(t_{f})\Delta
x(t_{f}),$

which can be rewritten as

$\displaystyle\big{(}I+\Psi(t_{f})\Delta\big{)}x(t_{f})=\mathrm{e}^{t_{f}\Lambda}x_{0}.$
(22)

Using the notation $H(t_{f})$ in (8), equation (22) is rewritten as

$\displaystyle x(t_{f})=H^{-1}(t_{f})\mathrm{e}^{t_{f}\Lambda}x_{0}.$ (23)

If the game has a unique open-loop Nash equilibrium, then (22) is satisfied
for any arbitrary $x_{0}$ and $x(t_{f})$. Equivalently, if matrix
$H^{-1}(t_{f})$ has an inverse for any arbitrary $x(t_{f})$, the unique
equilibrium actions exist and could be calculated for all $t\in[0,t_{f}]$. By
substituting (23) into (15) and re-substituting (15) in (13), I obtain (6).

Substituting $\lambda(t_{f})$ from (21) and then re-submitting $x(t_{f})$ from
(23), I have

$\displaystyle x(t)$
$\displaystyle=\mathrm{e}^{t\Lambda}x_{0}-\Psi(t)\lambda(t_{f})$
$\displaystyle=\mathrm{e}^{t\Lambda}x_{0}-\Psi(t)\Delta x(t_{f})$ (24)
$\displaystyle=\mathrm{e}^{t\Lambda}x_{0}-\Psi(t)\Delta
H^{-1}(t_{f})\mathrm{e}^{t_{f}\Lambda}x_{0},$

or in its final form (7). This concludes the proof. ∎

######  Remark 1.

From the definition of the Laplacian matrix $L_{i}$ in (4), it can be easily
figured out that all nonzero elements of $L_{i}$ are only in the $i$th row and
column. Moreover, the nonzero elements of the $i$th row and column are at the
indices $j\in\mathcal{N}_{i}$. Based on this structure, it can be deduced that
the matrix product $L_{i}H^{-1}(t_{f})\mathrm{e}^{t_{f}\Lambda}x_{0}$ only
requires the initial $x_{j0}$ entries for $\forall j\in\mathcal{N}_{i}$ in
$x_{0}$. Thus, I can draw the conclusion that the equilibrium actions (6) are
distributed in the sense that each agent only uses local information of her
own and her graph neighbors without using any global information of the
communication graph.

###  III-A Global Optimal

In the game-based optimization (5), each individual by minimizing her cost
function attains the locally optimal Nash equilibrium. However, the game
equilibrium in general does not correspond to the global social optimum, which
minimizes the sum of all costs.

Let $R=\mathrm{diag}(r_{1},\cdots,r_{n})$,
$B=\mathrm{diag}(b_{1},\cdots,b_{n})$, and
$u(t)=[u_{1}(t),\cdots,u_{n}(t)]^{\top}$. The global optimization problem for
the social network is

$\displaystyle\min_{u}~{}J(u(t))=x^{\top}(t_{f})Lx(t_{f})+\int_{0}^{t_{f}}u^{\top}(t)Ru(t)~{}\mathrm{dt},$
(25) $\displaystyle{s.t.}$ $\displaystyle\dot{x}(t)=\Lambda x(t)+Bu(t),\quad
x(0)=x_{0}.$

Define the Hamiltonian

$\mathcal{H}=u^{\top}(t)Ru(t)+\lambda^{\top}(t)\Big{(}\Lambda
x(t)+Bu(t)\Big{)}.$ (26)

Applying the necessary conditions on the Hamiltonian yield

$\displaystyle u(t)=-R^{-1}B^{\top}\lambda(t),$ (27)
$\displaystyle\dot{\lambda}(t)=-\Lambda^{\top}\lambda(t),\quad\lambda(t_{f})=Lx(t_{f}).$
(28)

The solution of (28) is

$\lambda(t)=\mathrm{e}^{(t_{f}-t)\Lambda^{\top}}\lambda(t_{f}).$ (29)

Substituting this solution in (27) yields the optimal control actions

$u(t)=-R^{-1}B^{\top}\mathrm{e}^{(t_{f}-t)\Lambda^{\top}}\lambda(t_{f}).$ (30)

Substituting (27) into the opinion dynamics in (25) and then using (27), I get

$\displaystyle\dot{x}(t)$ $\displaystyle=\Lambda
x(t)-BR^{-1}B^{\top}\lambda(t)$ $\displaystyle=\Lambda
x(t)-BR^{-1}B^{\top}\mathrm{e}^{(t_{f}-t)\Lambda^{\top}}\lambda(t_{f}).$ (31)

The solution of (III-A) at $t$ is given by

$x(t)=\mathrm{e}^{t\Lambda}x_{0}-\hat{\Psi}(t)\lambda(t_{f}),$ (32)

where

$\hat{\Psi}(t)=\int_{0}^{t}\mathrm{e}^{(t-\tau)\Lambda}BR^{-1}B^{\top}\mathrm{e}^{(t-\tau)\Lambda^{\top}}~{}\mathrm{d\tau}.$

Substituting $\lambda(t_{f})$ from (28) then into (32) at $t_{f}$ yields

$\displaystyle
x(t_{f})=\mathrm{e}^{t_{f}\Lambda}x_{0}-\hat{\Psi}(t_{f})Lx(t_{f}),$

which can be rewritten as

$\displaystyle\big{(}I+\hat{\Psi}(t_{f})L\big{)}x(t_{f})=\mathrm{e}^{t_{f}\Lambda}x_{0},$

or equivalently,

$\displaystyle x(t_{f})=\hat{H}^{-1}(t_{f})\mathrm{e}^{t_{f}\Lambda}x_{0},$
(33)

where

$\hat{H}(t_{f})=I+\hat{\Psi}(t_{f})L.$

Finally, the associated opinion trajectory with optimal control actions for
the global opinion formation problem in (25) is given by

$\displaystyle x(t)$
$\displaystyle=\big{(}\mathrm{e}^{t\Lambda}-\hat{\Psi}(t)L\hat{H}^{-1}(t_{f})\mathrm{e}^{t_{f}\Lambda}\big{)}x_{0}.$
(34)

Substituting $\lambda(t_{f})$ from (28) into (30) and then using (33), the
control actions are

$\displaystyle u(t$
$\displaystyle)=-R^{-1}B^{\top}\mathrm{e}^{(t_{f}-t)\Lambda^{\top}}L\hat{H}^{-1}(t_{f})\mathrm{e}^{t_{f}\Lambda}x_{0}.$
(35)

The nonsingularity of $\hat{H}(t_{f})$ is a prerequisite for (34) and (35).

##  IV Simulation Results

In this section, I apply the theoretical results to the well-known social
network of Zachary’s Karate Club [18]. With its underlying social graph given
in Fig. 1, Zachary’s Karate Club network, created by Wayne Zachary in 1977, is
a social network between the members of a karate club at a US university. This
network has 34 nodes and 78 edges, where the nodes and edges represent the
club members and their mutual friendships. I show the consequences of the
execution of game strategies in Theorem 1 on the evolution of opinions for
Zachary’s Karate Club network.

Figure 1: Zachary’s Karate Club network.

The simulations are carried out for the time interval of $t_{f}=10$ and the
initial opinion vector $x_{0}$ for each experiment is selected from a uniform
distribution in $(-1.5,-0.5)\cup(0.5,1.5)$. Such a distribution of the initial
opinions shows that the population’s initial opinions are split into two
groups or clusters. The evolution of club members’ opinions under the
continuous time HK model (2) prior to any optimization (i.e., $u_{i}(t)=0$) is
according to $x(t)=\mathrm{e}^{t\Lambda}x_{0}$. The corresponding opinion
trajectories are illustrated in Fig. 2. As it is seen, the club members’
opinions reach nearly a consensus about the average opinion in the network. I
should note that a consensus is realized at a long enough horizon length.

Figure 2: Evolution of opinion trajectories in the HK model for Zachary’s
network ($u_{i}(t)=0~{}\forall i\in\mathcal{V}$).

For optimizing the opinions, in the continuous-time HK model (2), I let
$b_{i}=1$. For Zachary’s Karate Club social network, the opinion trajectories
as a result of optimization (5) and (25) are shown in Fig. 3, by the green and
red colors, respectively. According to this figure, the non-stubborn network
under the framework of global optimization reaches a consensus on the average
opinions for $r_{i}=1$. Therefore, a consensus is the global social optimum
norm. Using their game strategies, the network members’ final opinions have
moved toward each other, but they have not reached a consensus. By imposing
large $r_{i}$, the network members seek to minimize their influence effort
rather than their disagreement with others. This is because, given a
sufficiently large $r_{i}$, the related optimization is reduced to the
minimization of the influence effort term. As it is seen from Fig. 3, for a
relatively large $r_{i}=20~{}\forall i\in\mathbf{V}$, the final opinion
trajectories have come closer to each other, but not as much as in the case
for $r_{i}=1$.

Figure 3: Opinion trajectories associated with the game strategies and global
optimal actions for the non-stubborn Zachary’s network, represented by the
green and red colors, respectively.

The finite horizon length is an important factor in the formation of final
opinions. From Fig. 1 and Fig. 3, it is seen that the optimization reduces the
finite horizon length for final opinion formation significantly. Fig. 4 shows
that a consensus emerges from the optimization of opinions at a long horizon
length.

Figure 4: A consensus of final opinions is on the horizon.

##  V Conclusions

This paper studied optimizing opinions of the Hegselmann-Krause model in a
game-theoretic framework. The execution of game equilibrium strategies in the
well-known Zachary’s Karate Club social network with non-stubborn individuals
showed that the opinions came close to each other but that a consensus of
final opinions did not emerge. However, the club members might be stubborn
individuals whose prejudices influence the formation of their final opinions.
A future research direction could be optimizing opinions, considering the
stubbornness of individuals in the game-theoretic analysis of social networks.
Feedback Nash equilibrium strategies for time-varying communication graphs and
confidence bound values can be investigated.

## Acknowledgment

This work was supported by SGS, VŠB - Technical University of Ostrava, Czech
Republic, under grant No. SP2023/12 “Parallel processing of Big Data X”.

## References

  * [1] H. M. Liang, C.C. Li, Y. C. Dong, Y. P. Jiang, “The fusion process of interval opinions based on the dynamic bounded confidence,” Information Fusion, vol. 29, 2016, pp. 112-119.
  * [2] N. E. Friedkin, “The problem of social control and coordination of complex systems in sociology: a look at the community cleavage problem,” IEEE Control Systems Magazine, vol. 35, no. 3, 2015, pp. 40-51.
  * [3] J. C. Dittmer, “Consensus formation under bounded confidence,” Nonlinear Analysis, vol. 47, no. 7, 2001, pp. 4615–4621.
  * [4] R. Hegselmann, U. Krause, “Opinion dynamics and bounded confidence models, analysis, and simulation,” Journal of Artificial Societies and Social Simulation, vol. 5, no. 3, 2002, pp. 1–33.
  * [5] G. Chen, W. Su, S. Ding, Y. Hong, “Heterogeneous Hegselmann–Krause dynamics with environment and communication noise,” IEEE Transactions on Automatic Control, vol. 65, no. 8, 2020, pp. 3409-3424.
  * [6] J. Lorenz, “Heterogeneous bounds of confidence: meet, discuss and find consensus!,” Complexity, vol. 15, no. 4, 2010, pp. 43–52.
  * [7] S.R. Etesami, T. Başar, “Game-theoretic analysis of the Hegselmann-Krause model for opinion dynamics in finite dimensions,” IEEE Transactions on Automatic Control, vol. 60, no. 7, 2015, pp. 1886-1897.
  * [8] E. Ghezelbash, M.J. Yazdanpanah, M. Asadpour, A. Yaghmaei, “Optimal selection of informed agents for influence opposition,” IEEE Transactions on Computational Social Systems, vol. 8, no. 1, 2021, pp. 20-31.
  * [9] D.S. Hunter , T. Zaman, “Optimizing opinions with stubborn agents,” Operations Research, vol. 70, no. 4, 2022, pp. 2119-2137.
  * [10] M.U.B. Niazi, A.B. Özgüler, “A differential game model of opinion dynamics: accord and discord as Nash equilibria,” Dynamic Games and Applications, vol. 11, 2021, pp. 137–160.
  * [11] R.A. Banez, H. Gao, L. Li, C. Yang, Z. Han, H. V. Poor, “Modeling and analysis of opinion dynamics in social networks using multiple-population mean field games,” IEEE Transactions on Signal and Information Processing over Networks, vol. 8, 2022, pp. 301-316.
  * [12] Z. Li, X. Chen, H. Yang, A. Szolnoki, “Game-theoretical approach for opinion dynamics on social networks,” Chaos, vol. 32, 2022, pp. 073117.
  * [13] A. Yıldız, A. B. Özgüler, “Opinion dynamics of stubborn agents under the presence of a troll as differential game,” Turkish Journal of Electrical Engineering & Computer Sciences, vol. 29, no. 7, 2021, pp. 23.
  * [14] S. Babakhanbak, K. Kavousi, F. Farokhi, “Application of a time delay neural network for predicting positive and negative links in social networks,” Turkish Journal of Electrical Engineering & Computer Sciences, vol. 24, no. 4, 2016, pp. 56.
  * [15] J. Engwerda, LQ Dynamic Optimization and Differential Games. John Wiley & Sons, 2005.
  * [16] E. V. Damme, Refinements of the Nash Equilibrium Concept, Lecture Notes in Economics and Mathematical Systems, Springer Berlin Heidelberg, 2012.
  * [17] H. B. Jond, V. V. Nabiyev, “On the finite horizon Nash equilibrium solution in the differential game approach to formation control,” Journal of Systems Engineering and Electronics, vol. 30, no. 6, 2019, pp. 1233-1242.
  * [18] W. W. Zachary, “An information flow model for conflict and fission in small groups,” Journal of Anthropological Research, vol. 33, no. 4, 1977, pp. 452-473.


# Azimuthal equatorial flows in spherical coordinates with discontinuous
stratification

Calin I. Martin

###### Abstract

We are concerned here with an exact solution to the governing equations for
geophysical fluid dynamics in spherical coordinates which incorporates
discontinuous fluid stratification. This solution represents a steady,
purely–azimuthal equatorial flow with an associated free-surface and an
interface separating two fluid regions, each of which having its own
continuous distribution of density. However, the two density functions do not
match along the interface. Following the derivation of the solution we
demonstrate that there is a well-defined relationship between the imposed
pressure at the free-surface and the resulting distortion of the surface’s
shape. Moreover, imposing the continuity of the pressure along the interface
generates an equation that describes (implicitly) the shape of the interface.
Interestingly, it turns out that the interface defining function has infinite
regularity.

Mathematics Subject Classification: 35Q31, 35Q35, 35Q86, 35R35, 76E20

Keywords: Azimuthal flows, discontinuous density, spherical coordinates,
Coriolis force, implicit function theorem.

##  1 Introduction

The Pacific Ocean displays within a band of 150 km on each side of the Equator
and extending longitudinally over about 16.000 km some outstanding features.
First and foremost we would like to point out a pronounced density
stratification (greater than anywhere else in the ocean [20]) whose hallmark
is the presence of a rather sharp interface (called pycnocline or thermocline)
that separates a shallow near-surface layer of relatively warm water from a
deep layer of colder and denser water. In fact, flow stratification (which in
our study varies with depth) represents an intrinsic trait of geophysical
water flows that particularly applies to large-scale ocean movements, cf. [21,
40, 50]. Greatly impacted by changes in temperature and salinity, density
fluctuations arise frequently in the ocean and lead to a layering of the flow:
fluid layers of different densities organize themselves so that the higher
densities are found below lower densities, cf. [22, 8, 13, 20, 38, 47]. This
vertical layering ushers an expected bundle of properties which considerably
influences the velocity field. As a result internal gravity waves may appear,
cf. Theorem 3.3.

On the other hand, the Coriolis forces (created by the Earth’s rotation)
together with the prevailing westward wind pattern generate an underlying
current field that displays flow-reversal [3, 16]: the current field changes
(in a band of about $2^{\circ}$ latitude around the Pacific Equator) from a
westward flow near the surface to an eastward-flowing jet called the
Equatorial Undercurrent (EUC) whose core resides more or less on the
thermocline. An interesting aspect in the equatorial ocean dynamics is that
the change of sign of the Coriolis force across the Equator turns the latter
into an actual waveguide that enables azimuthal flow propagation which is
symmetric with respect to the Equator.

Comprehensive mathematical models that handle wave-current interactions and
internal waves in the presence of density stratification within the setting of
nonlinear geophysical governing equations are quite rare and of very recent
date [5, 8, 9, 10, 15, 17, 34]. The stratification issue alone, in spite of
being of utmost practical relevance, has remained unamenable to a rigorous
mathematical analysis until recently: for a selection of recent advances in
the field (in the scenario of two-dimensional gravity water waves without
Earth’s rotation effects) we refer the reader to [9, 10, 15, 18, 19, 25, 46,
26, 51, 52].

In this paper, we address the topic of density stratification of discontinuous
type for flows exhibiting vertical structure, internal waves and a preferred
propagation direction and so derive and, subsequently, analyze an exact and
explicit solution which describes in a rotating frame a stratified inviscid
and incompressible azimuthal equatorial water flow which presents a free
surface and an interface. We would like to point out that explicit and/or
exact solutions in fluid dynamics represent a very rare occurrence. They are
important not only from a mathematical perspective but also for providing a
basis for more obtaining more elaborate and physically meaningful solutions by
employing perturbative or asymptotic methods [37].

Mindful of the previous aspects we aim at furthering a line of work put
forward by Constantin & Johnson concerning mathematical analyses of
geophysical water flows. First, Constantin & Johnson constructed in [11, 12]
in terms of spherical coordinates exact solutions to the full GFD governing
equations. These solutions represent purely-azimuthal, depth-varying flows,
can be chosen to model both the equatorial undercurrent (EUC), and the
Antarctic Circumpolar Current (ACC), respectively, and describe purely
homogeneous flows without stratification. Recently, Henry & Martin presented
equatorial flow solutions that allow for continuous stratification that varies
linearly with depth [27, 28] as well as of most general type [29, 30]. Martin
& Quirchmayr [41] devised solutions flows modelling the ACC exhibiting a
general fluid stratification that varies both with depth and latitude. The
special trait of the solution that we derive consists in its ability to
accommodate a discontinuous density stratification that varies with depth.
More precisely, we allow a vertical layering as described in the beginning,
with two layers of different, non-constant densities, where the denser layer
sits below the less dense one (stable stratification): thus an interface
arises which resembles the role played by the thermocline in the setting of
observed equatorial flows.

Subsequent to the presentation of the equations of motion in Section 2 we
derive in Section 3 explicit formulas for the velocity field and for the
pressure function. We then exploit the dynamic boundary condition and find a
relation between the imposed pressure on the surface and the resulting surface
distortion. Moreover, from the balance of forces at the interface between the
two layers, we derive another equation for the interface defining function.
Although the two relations which define implicitly the free surface and the
interface, respectively, appear in a convoluted form, they have the redeeming
feature of being amenable to a functional analytic study by means of the
implicit function theorem. We present in Section 4 results which show that our
exact solution manifests acceptable properties: a growth in pressure along the
surface leads to a decline in height of the free surface. Moreover, we also
show that the interface between the two layers has infinite regularity.

##  2 Equations of motion

The governing equations for geophysical fluid dynamics (GFD) are formulated in
a spherical coordinate system which is fixed at a point on the Earth’s surface
as follows. We work in spherical coordinates, denoted $(r,\theta,\varphi)$,
where $r$ is the distance from the centre of the earth, $\theta$ (with
$0\leq\theta\leq\pi$) is the polar angle, and $\varphi$ (with
$0\leq\varphi<2\pi$) is the azimuthal angle, see Figure 1. The location of the
North and South poles are at $\theta=0,\pi$, respectively, while the Equator
is situated at $\theta=\frac{\pi}{2}$. The unit vectors in this
$(r,\theta,\varphi)$ system are
$(\mathbf{e}_{r},\mathbf{e}_{\theta},\mathbf{e}_{\varphi})$, respectively, and
the corresponding velocity components are $({\bf w},{\bf v},{\bf u})$;
$\mathbf{e}_{\varphi}$ points from West to East, and $\mathbf{e}_{\theta}$
points from North to South. A noteworthy aspect is that the fluid domain is
stratified where the stratification is brought about by the changes in the
(discontinuous) density. More precisely, denoting with $R\approx 6378$ km the
Earth’s radius and with $r_{0}>0,r_{1}$ some constants, we assume that the
flow consists of a top layer

$D_{top}:=\\{(r,\theta,\varphi):R+r_{1}+h(\theta,\varphi)\leq r\leq
R+r_{0}+k(\theta,\varphi)\\},$

that lies above a bottom layer

$D_{bottom}:=\\{(r,\theta,\varphi):R+d(\theta,\varphi)\leq r\leq
R+r_{1}+h(\theta,\varphi)\\}.$

While above $d$ is a given function, $h$ and $k$ are unknowns of the problem
and denote the interface and the free surface defining functions,
respectively.

Figure 1: The spherical coordinate system as presented in Section 2. The
variable $r$ represents the distance from the origin, $\theta$ is the polar
angle ($\pi/2-\theta$ being the angle of latitude), and $\varphi$ is the angle
of longitude (azimuthal angle).

To draw attention to the layering aspect of the fluid domain we will use the
index $1$ as a label for the upper layer. Whenever we refer to the overall
physical variable without specification of the layer, we shall use bold face
symbol. For example, the density $\bm{\rho}$ is assumed to be discontinuous
across $r=R+h(\theta,\varphi)$ and is distributed in the two layers as

$\bm{\rho}(r)=\left\\{\begin{array}[]{l}\rho_{1}(r)\,\,{\rm
for}\,\,(r,\theta,\varphi)\in D_{top},\\\ \rho(r)\,\,{\rm
for}\,\,(r,\theta,\varphi)\in D_{bottom},\end{array}\right.$ (1)

being a differentiable function of $r$ in each of $D_{top}$ and $D_{bottom}$.
Likewise, the components of the velocity field with respect to the orthonormal
system $\\{\mathbf{e}_{r},\mathbf{e}_{\theta},\mathbf{e}_{\varphi}\\}$ are
denoted with $({\bf w},{\bf v},{\bf u})$. That is,

$({\bf w},{\bf v},{\bf
u})=\left\\{\begin{array}[]{l}(w_{1},v_{1},u_{1})\,\,{\rm in}\,\,D_{top},\\\
(w,v,u)\,\,{\rm in}\,\,D_{bottom}.\end{array}\right.$ (2)

We present now the equations of motion in a coordinate system with its origin
at the centre of the sphere, where the effects of the earth’s rotation are
incorporated into our system. In order to do this we associate
$\\{\mathbf{e}_{r},\mathbf{e}_{\theta},\mathbf{e}_{\varphi}\\}$ with a point
fixed on the sphere which is rotating about its polar axis. Consequently, the
governing equations will contain Coriolis $2{\bm{\Omega}}\times\bf{U}$ terms
and centripetal ${\bm{\Omega}}\times({\bm{\Omega}}\times\bf{r})$ terms, with
${\bf r}=r\mathbf{e}_{r}$, ${\bf U}={\bf w}\mathbf{e}_{r}+{\bf
v}\mathbf{e}_{\theta}+{\bf u}\mathbf{e}_{\varphi}$, and
${\bm{\Omega}}=\varOmega(\mathbf{e}_{r}\cos\theta-\mathbf{e}_{\theta}\sin\theta),$
where $\varOmega\approx 7.29\times 10^{-5}$ rad s-1 is the constant rate of
rotation of the earth. The contribution of the Coriolis and centripetal
acceleration is therefore given by

$2\varOmega(-{\bf u}\sin\theta\mathbf{e}_{r}-{\bf
u}\cos\theta\mathbf{e}_{\theta}+({\bf w}\sin\theta+{\bf
v}\cos\theta)\mathbf{e}_{\varphi})-r\varOmega^{2}(\sin^{2}\theta\mathbf{e}_{r}+\sin\theta\cos\theta\mathbf{e}_{\theta}).$

Therefore, the GFD governing equations, in a spherical coordinate system with
its origin at the centre of the sphere, are given by the Euler equation

$\begin{split}{\bf w}_{t}+{\bf w}{\bf w}_{r}+\frac{{\bf v}}{r}{\bf
w}_{\theta}+\frac{{\bf u}}{r\sin\theta}{\bf w}_{\varphi}-\frac{1}{r}({\bf
v}^{2}+{\bf u}^{2})&-2\varOmega{\bf
u}\sin\theta-r\varOmega^{2}\sin^{2}\theta\\\ &=-\frac{{\bf
p}_{r}}{\bm{\rho}}+\mathfrak{F}_{r}\\\ \mathbf{v}_{t}+{\bf w}{\bf
v}_{r}+\frac{{\bf v}}{r}{\bf v}_{\theta}+\frac{{\bf u}}{r\sin\theta}{\bf
v}_{\varphi}+\frac{1}{r}({\bf w}{\bf v}-{\bf u}^{2}\cot\theta)&-2\varOmega{\bf
u}\cos\theta-r\varOmega^{2}\sin\theta\cos\theta\\\ &=-\frac{{\bf
p}_{\theta}}{\bm{\rho}r}+\mathfrak{F}_{\theta}\\\ {\bf u}_{t}+{\bf w}{\bf
u}_{r}+\frac{{\bf v}}{r}{\bf u}_{\theta}+\frac{{\bf u}}{r\sin\theta}{\bf
u}_{\varphi}+\frac{1}{r}({\bf w}{\bf u}+{\bf v}{\bf
u}\cot\theta)&+2\varOmega{\bf w}\sin\theta+2\varOmega{\bf v}\cos\theta\\\
&=-\frac{{\bf
p}_{\varphi}}{\bm{\rho}r\sin\theta}+\mathfrak{F}_{\varphi},\end{split}$ (3)

and the equation of mass conservation

$\frac{1}{r^{2}}\frac{\partial}{\partial r}(r^{2}\bm{\rho}{\bf
w})+\frac{1}{r\sin\theta}\frac{\partial}{\partial\theta}(\bm{\rho}{\bf
v}\sin\theta)+\frac{1}{r\sin\theta}\frac{\partial}{\partial\varphi}(\bm{\rho}{\bf
u})=0.$ (4)

Here ${\bf p}(r,\theta,\varphi)$ denotes the pressure in the fluid,
$\mathfrak{F}=\mathfrak{F}_{r}\mathbf{e}_{r}+\mathfrak{F}_{\theta}\mathbf{e}_{\theta}+\mathfrak{F}_{\varphi}\mathbf{e}_{\varphi}$
is the body-force vector. The governing equations are supplemented by
associated boundary conditions as follows. At the free-surface we require the
dynamic condition involving the surface pressure

$p_{1}=P_{1}(\theta,\varphi),$ (5a) and the kinematic condition
$w_{1}=\frac{v_{1}}{r}\frac{\partial
k}{\partial\theta}+\frac{u_{1}}{r\sin\theta}\frac{\partial
k}{\partial\varphi},$ (5b) to hold.

At the interface, $r=R+h(\theta,\varphi)$, we ask that the normal components
of the velocity fields from the upper and lower layer, respectively, are the
same. The latter is equivalent with the condition

$\begin{split}(w\mathbf{e}_{r}+v\mathbf{e}_{\theta}+u\mathbf{e}_{\varphi})&\cdot\left(\mathbf{e}_{r}-\frac{h_{\theta}}{r}\mathbf{e}_{\theta}-\frac{h_{\varphi}}{r\sin\theta}\mathbf{e}_{\varphi}\right)\\\
&=(w_{1}\mathbf{e}_{r}+v_{1}\mathbf{e}_{\theta}+u_{1}\mathbf{e}_{\varphi})\cdot\left(\mathbf{e}_{r}-\frac{h_{\theta}}{r}\mathbf{e}_{\theta}-\frac{h_{\varphi}}{r\sin\theta}\mathbf{e}_{\varphi}\right).\end{split}$
(5c)

Moreover, to ensure the balance of forces at the interface, we also require
that the pressure from the upper layer coincides with the pressure from the
bottom layer, that is,

$p(R+h(\theta,\varphi),\theta,\varphi)=p_{1}(R+h(\theta,\varphi),\theta,\varphi).$
(5d)

At the bottom of the ocean, which is an impermeable, solid boundary described
by the equation $r=d(\theta,\varphi)$, the associated kinematic condition is

$w=\frac{v}{r}\frac{\partial
d}{\partial\theta}+\frac{u}{r\sin\theta}\frac{\partial d}{\partial\varphi}.$
(5e)

##  3 Explicit and exact solutions

We set out to derive exact solutions to the equations of motion and boundary
conditions from Section 2. In fact, the velocity field and the pressure
function are given explicitly. However, due to intricacies like the usage of
spherical coordinates, the complexity of stratification and the very involved
formulas for the pressure and velocity, explicite formulas for the free
surface and the interface remain elusive. An alleviation of this aspect is
represented by a Bernoulli-type relation between the imposed pressure at the
surface and the resulting surface distorsion; an implicit formula for the
interface defining function being given by the balance of forces at the
interface (5d). For a selective list of papers presenting exact solutions
pertaining to geophysical fluid dynamics we refer the reader to [1, 4, 6, 7,
11, 12, 13, 14, 16, 23, 24, 31, 32, 33, 39, 42, 43, 44, 45, 48].

Our objective is to find solutions to equations (3), (4) and (5) that
represent purely-azimuthal steady flows with no variation in the azimuthal
direction. Hence, the velocity field satisfies $w=v=w_{1}=v_{1}=0$ and
$u=u(r,\theta),\,u_{1}=u_{1}(r,\theta)$. Moreover, the other unknows of the
problem are characterized by the properties
$p=p(r,\theta),\,p_{1}=p_{1}(r,\theta),\,h=h(\theta),\,k=k(\theta)$. We notice
that a flow with the previous features automatically satisfies the boundary
and interface conditions (5) as well as the equation of mass conservation (4).
As for the Euler equations, they become equivalent with

$\left\\{\begin{array}[]{rcl}-\frac{{\bf u}^{2}}{r}-2\varOmega{\bf
u}\sin\theta-r\varOmega^{2}\sin^{2}\theta&=&-\frac{1}{\rho}p_{r}-g\\\ &&\\\
-\frac{{\bf u}^{2}}{r}\cot\theta-2\varOmega{\bf
u}\cos\theta-r\varOmega^{2}\sin\theta\cos\theta&=&-\frac{1}{\rho
r}p_{\theta}\\\ &&\\\
0&=&-\frac{1}{\rho}\frac{1}{r\sin\theta}p_{\varphi}\end{array},\right.$ (6)

where we assumed that the external body-force is due to gravity alone, giving
the body-force vector as $-g\mathbf{e}_{r}$. Eliminating the pressure from (6)
we obtain the equation

$\partial_{\theta}\left(\frac{\rho(r)({\bf u}(r,\theta)+\Omega
r\sin\theta)^{2}}{r}\right)-\partial_{r}\left(\rho(r)({\bf u}(r,\theta)+\Omega
r\sin\theta)^{2}\cot\theta\right)=0.$ (7)

To solve (7) we appeal to the method of characteristics. An analysis similar
to the one in [29] gives for the velocity field the formula

${\bf u}(r,\theta)=\left\\{\begin{aligned} u(r,\theta)=-\Omega
r\sin\theta+\frac{F(r\sin\theta)}{\sqrt{\rho(r)}}\quad{\rm if}\quad
R+d(\theta)\leq r\leq R_{1}+h(\theta),\\\ u_{1}(r,\theta)=-\Omega
r\sin\theta+\frac{F_{1}(r\sin\theta)}{\sqrt{\rho_{1}(r)}}\quad{\rm if}\quad
R_{1}+h(\theta)\leq r\leq R_{0}+k(\theta),\end{aligned}\right.$ (8)

where $R_{0}:=R+r_{0},R_{1}:=R+r_{1}$, and $x\rightarrow F(x),\,x\rightarrow
F_{1}(x)$ are arbitrary real-valued functions and $r_{0},r_{1}$ are some
positive constants. Further, inserting the expressions of $w$ and $w_{1}$ from
(8) in (6) we have

${\bf p}_{r}=\left\\{\begin{aligned}
p_{r}(r,\theta)=\frac{F^{2}(r\sin\theta)}{r}-g\rho(r)\quad{\rm if}\quad
R+d(\theta)\leq r\leq R_{1}+h(\theta),\\\
p_{1,r}(r,\theta)=\frac{F_{1}^{2}(r\sin\theta)}{r}-g\rho_{1}(r)\quad{\rm
if}\quad R_{1}+h(\theta)\leq r\leq R_{0}+k(\theta),\end{aligned}\right.$ (9)

and

${\bf p}_{\theta}=\left\\{\begin{aligned}
p_{\theta}(r,\theta)=F^{2}(r\sin\theta)\cot\theta\quad{\rm if}\quad
R+d(\theta)\leq r\leq R_{1}+h(\theta),\\\
p_{1,\theta}(r,\theta)=F_{1}^{2}(r\sin\theta)\cot\theta\quad{\rm if}\quad
R_{1}+h(\theta)\leq r\leq R_{0}+k(\theta).\end{aligned}\right.$ (10)

From integration with respect to $r$ we obtain for all
$r\in[R+d(\theta),R_{1}+h(\theta)]$

$p(r,\theta)=\int_{(R+d(\theta))\sin\theta}^{r\sin\theta}\frac{F^{2}(y)}{y}dy-g\int_{R+d(\theta)}^{r}\rho(\tilde{r})d\tilde{r}+f(\theta),$
(11)

where

$f(\theta)=\int_{R+d(\pi/2)}^{(R+d(\theta))\sin\theta}\frac{F^{2}(y)}{y}dy-g\int_{R+d(\pi/2)}^{R+d(\theta)}\rho(y)dy.$
(12)

For the upper layer we have

$p_{1}(r,\theta)=\int_{(R_{1}+h(\theta))\sin\theta}^{r\sin\theta}\frac{F_{1}^{2}(y)}{y}dy-g\int_{R_{1}+h(\theta)}^{r}\rho_{1}(\tilde{r})d\tilde{r}+c_{1}(h,\theta)$
(13)

for all $r\in[R_{1}+h(\theta),R_{0}+k(\theta)]$, where

$\displaystyle c_{1}(h,\theta)=$
$\displaystyle\int_{\pi/2}^{\theta}F_{1}^{2}\big{(}(R_{1}+h(\tilde{\theta}))\sin\tilde{\theta}\big{)}\left[\cot\tilde{\theta}+\frac{h^{\prime}(\tilde{\theta})}{R_{1}+h(\tilde{\theta})}\right]\,d\tilde{\theta}$
(14)
$\displaystyle-g\int_{\pi/2}^{\theta}\rho_{1}(R_{1}+h(\tilde{\theta}))h^{\prime}(\tilde{\theta})\,d\tilde{\theta}.$

We pass now to the (implicit) determination of the two interfaces: the free
surface and the interface separating the two layers $D_{bottom}$ and
$D_{top}$. From the dynamic condition on the surface (5a) we obtain

$P_{1}(\theta)=\int_{(R_{1}+h(\theta))\sin\theta}^{(R_{0}+k(\theta))\sin\theta}\frac{F_{1}^{2}(y)}{y}dy-g\int_{R_{1}+h(\theta)}^{R_{0}+k(\theta)}\rho_{1}(\tilde{r})d\tilde{r}+c_{1}(h,\theta).$
(15)

Setting $h=0$ and $k=0$ in the previous formula we obtain the quantity

$P^{0}_{1}(\theta)=P_{1}(h\equiv
0)=\int_{R_{1}\sin\theta}^{R_{0}\sin\theta}\frac{F_{1}^{2}(y)}{y}dy-g\int_{R_{1}}^{R_{0}}\rho_{1}(\tilde{r})d\tilde{r}+c_{1}(0,\theta),$
(16)

which stands for the pressure required to maintain an undisturbed free surface
and an undisturbed interface following the curvature of the Earth.
The balance of forces at the interface, encoded in the relation

$p(R_{1}+h(\theta),\theta)=p_{1}(R_{1}+h(\theta),\theta),$ (17)

will allow us to implicitly determine the shape of the interface. Equation
(17) will be studied in a functional analytic framework by means of
nondimensionalization. More precisely, setting

$\mathcal{h}(\theta):=\frac{h(\theta)}{R_{1}}$ (18)

we see that (17) can be written as $\mathcal{G}(\mathcal{h})=0$, where

$\displaystyle\mathcal{G}(\mathcal{h})(\theta)=$ (19)
$\displaystyle\frac{1}{P_{atm}}\left(\int_{(R+d(\theta))\sin\theta}^{(1+\mathcal{h}(\theta))R_{1}\sin\theta}\frac{F^{2}(y)}{y}dy-g\int_{R+d(\theta)}^{R_{1}(1+\mathcal{h}(\theta))}\rho(\tilde{r})d\tilde{r}+f(\theta)\right)$
$\displaystyle-\frac{1}{P_{atm}}\int_{\pi/2}^{\theta}F_{1}^{2}\big{(}(1+\mathcal{h}(\tilde{\theta}))R_{1}\sin\tilde{\theta}\big{)}\left[\cot\tilde{\theta}+\frac{\mathcal{h}^{\prime}(\tilde{\theta})}{1+\mathcal{h}(\tilde{\theta})}\right]\,d\tilde{\theta}$
$\displaystyle+\frac{gR_{1}}{P_{atm}}\int_{\pi/2}^{\theta}\rho_{1}\big{(}R_{1}(1+\mathcal{h}(\tilde{\theta}))\big{)}\mathcal{h}^{\prime}(\tilde{\theta})\,d\tilde{\theta},$

where $f(\theta)$ is given in (12). Also, the implicit determination of the
free surface will be devised from the dynamic condition (5a), the latter being
written in the abstract form

$\mathcal{F}(\mathcal{h},\mathcal{k},\mathcal{P}_{1})=0,$ (20)

with $\mathcal{F}$ acting from the Banach space
$C^{1}\big{(}\left(\left[\pi/2,\pi/2+\varepsilon\right]\right)\times
C\left(\left[\pi/2,\pi/2+\varepsilon\right]\right)\times
C\left(\left[\pi/2,\pi/2+\varepsilon\right]\right)\big{)}$ into itself, where
$\varepsilon=0.016$ is suitable for flows in the equatorial region and
demarcates a strip of about 100 km about the Equator, cf. [11]. Above
$\mathcal{k}(\theta):=\frac{k(\theta)}{R_{0}},\,\,\mathcal{P}_{1}(\theta):=\frac{P(\theta)}{P_{atm}}$
and

$\begin{split}&\mathcal{F}(\mathcal{h},\mathcal{k},\mathcal{P}_{1})\\\
&=\frac{1}{P_{atm}}\left(\int_{(1+h(\theta))R_{1}\sin\theta}^{(1+\mathcal{k}(\theta))R_{0}\sin\theta}\frac{F_{1}^{2}(y)}{y}dy-g\int_{R_{1}(1+\mathcal{h}(\theta))}^{R_{0}(1+\mathcal{k}(\theta))}\rho_{1}(\tilde{r})d\tilde{r}+c_{1}(h,\theta)\right)-\mathcal{P}_{1}(\theta),\end{split}$
(21)

where $c_{1}$ is given by (14). We can at once identify one trivial solution
of (20), namely the one associated with the flow having the free surface and
the interface in an undisturbed state, following the curvature of the Earth.
That is, we have

$\mathcal{F}(0,0,\mathcal{P}_{1}^{0})=0,$

where

$\mathcal{P}_{1}^{0}=\frac{P_{1}^{0}}{P_{atm}}=\frac{1}{P_{atm}}\left(\int_{R_{1}\sin\theta}^{R_{0}\sin\theta}\frac{F_{1}^{2}(y)}{y}dy-g\int_{R_{1}}^{R_{0}}\rho_{1}(\tilde{r})\,d\tilde{r}+\int_{\pi/2}^{\theta}F_{1}^{2}(R_{1}\sin\tilde{\theta})\cot\tilde{\theta}\,d\tilde{\theta}\right).$

We recast now the problem of finding $h$ and $k$ as the problem of finding
solutions to the equation

$(\mathcal{G}(h),\mathcal{F}(\mathcal{h},\mathcal{k},P_{1}))=0.$ (22)

The previous task will be carried out by resorting to the implicit function
theorem. Therefore, we will need to compute the derivative of the map
$(\mathcal{G},\mathcal{F})$ with respect to $(h,k)$. To this end we compute
first

$\displaystyle\frac{c_{1}(s\mathcal{h},\theta)-c_{1}(0,\theta)}{s}=$
$\displaystyle\frac{1}{s}\int_{\pi/2}^{\theta}\left[F_{1}^{2}\big{(}(R_{1}+sR_{1}\mathcal{h}(\tilde{\theta}))\sin\tilde{\theta}\big{)}-F_{1}^{2}(R_{1}\sin\tilde{\theta})\right]\cot\tilde{\theta}\,d\tilde{\theta}$
$\displaystyle+\int_{\pi/2}^{\theta}\frac{\mathcal{h}^{\prime}(\tilde{\theta})}{1+s\mathcal{h}(\tilde{\theta})}F_{1}^{2}\big{(}(R_{1}+sR_{1}\mathcal{h}(\tilde{\theta}))\sin\tilde{\theta}\big{)}\,d\tilde{\theta}$
$\displaystyle-
gR_{1}\int_{\pi/2}^{\theta}\rho_{1}(R_{1}+sR_{1}\mathcal{h}(\tilde{\theta}))\mathcal{h}^{\prime}(\tilde{\theta})\,d\tilde{\theta}.$

Utilizing the mean value theorem for integrals we obtain

$\displaystyle\lim_{s\to 0}\frac{c_{1}(sh,\theta)-c_{1}(0,\theta)}{s}=$
$\displaystyle\int_{\pi/2}^{\theta}\left[(F_{1}^{2})^{\prime}(R_{1}\sin\tilde{\theta})R_{1}\mathcal{h}(\tilde{\theta})\cos\tilde{\theta}+\mathcal{h}^{\prime}(\tilde{\theta})F_{1}^{2}\big{(}R_{1}\sin\tilde{\theta}\big{)}\right]\,d\tilde{\theta}$
(23) $\displaystyle-
gR_{1}\rho_{1}(R_{1})\big{(}\mathcal{h}(\theta)-\mathcal{h}(\pi/2)\big{)}$
$\displaystyle=$
$\displaystyle[F_{1}^{2}(R_{1}\sin\theta)-gR_{1}\rho_{1}(R_{1})]\mathcal{h}(\theta),$

where, in the last step, we took into account that $\mathcal{h}(\pi/2)=0$.
Using (23) we obtain

$\displaystyle P_{atm}\cdot(\mathcal{G}_{\mathcal{h}}(0)\mathcal{h})(\theta)=$
$\displaystyle\lim_{s\to 0}\frac{\mathcal{G}(s\mathcal{h})-\mathcal{G}(0)}{s}$
(24) $\displaystyle=$ $\displaystyle\lim_{s\to
0}\frac{1}{s}\left[\int_{R_{1}\sin\theta}^{(R_{1}+sR_{1}\mathcal{h}(\theta))\sin\theta}\frac{F^{2}(y)}{y}\,dy-g\int_{R_{1}}^{R_{1}+sR_{1}\mathcal{h}(\theta)}\rho(\tilde{r})\,d\tilde{r}\right]$
$\displaystyle-\lim_{s\to
0}\frac{1}{s}\int_{\pi/2}^{\theta}\left(F_{1}^{2}((R_{1}+sR_{1}\mathcal{h}(\tilde{\theta}))\sin\tilde{\theta})-F_{1}^{2}(R_{1}\sin\tilde{\theta})\right)\cot\tilde{\theta}\,d\tilde{\theta}$
$\displaystyle-\lim_{s\to
0}\int_{\pi/2}^{\theta}F_{1}^{2}((R_{1}+sR_{1}\mathcal{h}(\tilde{\theta}))\sin\tilde{\theta})\frac{\mathcal{h}^{\prime}(\tilde{\theta})}{1+s\mathcal{h}(\tilde{\theta})}\,d\tilde{\theta}$
$\displaystyle+gR_{1}\lim_{s\to
0}\int_{\pi/2}^{\theta}\rho_{1}(R_{1}+sR_{1}\mathcal{h}(\tilde{\theta}))\mathcal{h}^{\prime}(\tilde{\theta})\,d\tilde{\theta}.$

Appealing to the mean value theorem for integrals we obtain

$(\mathcal{G}_{h}(0)h)(\theta)=\big{[}F^{2}(R_{1}\sin\theta)-F_{1}^{2}(R_{1}\sin\theta)-gR_{1}(\rho(R_{1})-\rho_{1}(R_{1}))\big{]}\mathcal{h}(\theta).$
(25)

###### Remark 3.1

Field data suggest that the sizes of the densities $\rho$ and $\rho_{1}$
relate as $\rho=\rho_{1}(1+\sigma)$, where $\tau\rightarrow\sigma(\tau)$ is a
positive function satisfying $\sigma=\mathcal{O}(10^{-3})$. Moreover, guided
by physical reality, we also assume that
$\rho_{1}=\mathcal{O}(10^{3})\,kg\cdot m^{-3}$. Accounting for the previous
considerations we find that

$\displaystyle F^{2}(R_{1}\sin\theta)$ $\displaystyle-
F_{1}^{2}(R_{1}\sin\theta)=$ (26)
$\displaystyle=\rho(R_{1})(u(R_{1},\theta)+\Omega
R_{1}\sin\theta)^{2}-\rho_{1}(R_{1})(u_{1}(R_{1},\theta)+\Omega
R_{1}\sin\theta)^{2}$ $\displaystyle=\rho_{1}(R_{1})(u^{2}-u_{1}^{2}+2\Omega
R_{1}(u-u_{1})\sin\theta)+\rho_{1}\sigma(u+\Omega R_{1}\sin\theta)^{2}$
$\displaystyle<\rho_{1}u^{2}+2\Omega
R_{1}\rho_{1}u\sin\theta+\rho_{1}\sigma(u+\Omega R_{1}\sin\theta)^{2}.$

Allowing for values of $w,\rho_{1}$, that bear relevance on physical grounds,
we have from above that

$F^{2}(R_{1}\sin\theta)-F_{1}^{2}(R_{1}\sin\theta)<R_{1}\,kg\cdot m^{-1}\cdot
s^{-2},$

while

$gR(\rho-\rho_{1})=gR_{1}\rho_{1}\sigma>4.9R_{1}\,kg\cdot m^{-1}\cdot s^{-2}.$

The previous two inequalities show that there is a constant $\alpha<0$ such
that

$F^{2}(R_{1}\sin\theta)-F_{1}^{2}(R_{1}\sin\theta)-gR_{1}(\rho(R_{1})-\rho_{1}(R_{1}))\leq\alpha$

holds true for all $\theta\in[0,\pi]$. The latter finding guarantees that the
map

$\mathcal{h}\rightarrow\mathcal{G}_{\mathcal{h}}(0)\mathcal{h}$

is a linear homeomorphism from
$C^{1}\left(\left[\pi/2,\pi/2+\varepsilon\right]\right)$ to
$C^{1}\left(\left[\pi/2,\pi/2+\varepsilon\right]\right)$.

We further compute

$\displaystyle\mathcal{F}_{\mathcal{h}}(0,0,P_{1}^{0})\mathcal{h}=$
$\displaystyle\lim_{s\to
0}\frac{\mathcal{F}(s\mathcal{h},0,P_{1}^{0})-\mathcal{F}(0,0,P_{1}^{0})}{s}$
(27) $\displaystyle=$
$\displaystyle\frac{-F_{1}^{2}(R_{1}\sin\theta)\mathcal{h}+gR_{1}\rho_{1}(R_{1})\mathcal{h}+F_{1}^{2}(R_{1}\sin\theta)\mathcal{h}-gR_{1}\rho_{1}(R_{1})\mathcal{h}}{P_{atm}}$
$\displaystyle=$ $\displaystyle 0.$

Finally, utilizing (8), we obtain that for all $\theta\in[0,\pi]$ it holds

$\displaystyle(\mathcal{F}_{\mathcal{k}}(0,0,P_{1}^{0})\mathcal{k})(\theta)=$
$\displaystyle\lim_{s\to
0}\frac{\mathcal{F}(0,s\mathcal{k},P_{1}^{0})-\mathcal{F}(0,0,P_{1}^{0})}{s}$
(28) $\displaystyle=$
$\displaystyle\frac{F_{1}^{2}(R_{0}\sin\theta)\mathcal{k}(\theta)-gR_{0}\rho_{1}(R_{0})\mathcal{k}(\theta)}{P_{atm}}$
$\displaystyle=$
$\displaystyle\frac{\rho_{1}(R_{0})}{P_{atm}}\left[(u(R_{0},\theta)+\Omega
R_{0}\sin\theta)^{2}-gR_{0}\right]\mathcal{k}(\theta).$

###### Remark 3.2

Arguing as in Remark 3.1 we conclude that there is $\beta<0$ such that
$u(R_{0},\theta)+\Omega R_{0}\sin\theta)^{2}-gR_{0}\leq\beta$ for all
$\theta\in[0,\pi]$. Consequently, we have that the map

$\mathcal{k}\rightarrow\mathcal{F}_{\mathcal{k}}(0,0,P_{1}^{0})\mathcal{k}$

is a linear homeomorphism from
$C\left(\left[\pi/2,\pi/2+\varepsilon\right]\right)$ to
$C\left(\left[\pi/2,\pi/2+\varepsilon\right]\right)$.

Hence

$(\mathcal{G},\mathcal{F})_{\mathcal{h},\mathcal{k}}(0,0,P_{1}^{0})=\begin{pmatrix}\mathcal{G}_{\mathcal{h}}(0,0,P_{1}^{0})&\mathcal{G}_{\mathcal{k}}(0,0,P_{1}^{0})\\\
\mathcal{F}_{\mathcal{h}}(0,0,P_{1}^{0})&\mathcal{F}_{\mathcal{k}}(0,0,P_{1}^{0})\end{pmatrix}=\begin{pmatrix}\mathcal{G}_{\mathcal{h}}(0,0,P_{1}^{0})&0\\\
0&\mathcal{F}_{\mathcal{k}}(0,0,P_{1}^{0})\end{pmatrix}$

is a linear homeomorphism from
$C^{1}\left(\left[\pi/2,\pi/2+\varepsilon\right]\right)\times
C\left(\left[\pi/2,\pi/2+\varepsilon\right]\right)$ to
$C^{1}\left(\left[\pi/2,\pi/2+\varepsilon\right]\right)\times
C\left(\left[\pi/2,\pi/2+\varepsilon\right]\right)$ by Remark 3.1 and Remark
3.2. Hence, by the implicit function theorem, we have the following existence
result concerning the free surface and the interface.

###### Theorem 3.3

For any sufficiently small perturbation $\mathcal{P}_{1}$ of
$\mathcal{P}_{1}^{0}$ there exists a unique tuple
$(\mathcal{h},\mathcal{k})\in
C^{1}\left(\left[\pi/2,\pi/2+\varepsilon\right]\right)\times
C\left(\left[\pi/2,\pi/2+\varepsilon\right]\right)$ such that (22) holds true.

We would like to conclude this section by illustrating the usefulness of the
exact solutions we found. More precisely, we will show that observed salient
features of equatorial flows (like a pronounced density stratification, giving
rise to the pycnocline, as well as the presence of depth-dependent current
fields) can be successfully captured by the setting provided in (8).

###### Remark 3.4

Firstly, we note that it is possible to determine the azimuthal flow velocity
at latitude $\theta$ from its value at the Equator ($\theta=\pi/2$): indeed,
it is a consequence of formula (8) that

$u(r,\theta)=u(r,\pi/2)\sin\theta+\frac{1}{\sqrt{\rho(r)}}(F(r\sin\theta)-F(r)\sin\theta)\,\,\,{\rm
in}\,\,\,\Omega_{bottom},$ (29)

and

$u_{1}(r,\theta)=u_{1}(r,\pi/2)\sin\theta+\frac{1}{\sqrt{\rho_{1}(r)}}(F_{1}(r\sin\theta)-F_{1}(r)\sin\theta)\,\,\,{\rm
in}\,\,\,\Omega_{top}.$ (30)

Moreover, the functions $F$ and $F_{1}$ can be reconstituted from the
equatorial part of the flow as

$F(r)=\sqrt{\rho(r)}(u(r,\pi/2)+\Omega r),\quad
F_{1}(r)=\sqrt{\rho_{1}(r)}(u_{1}(r,\pi/2)+\Omega r).$ (31)

Next we make use of formulas (29),(30) and (31) to show that the explicit
solutions derived here are relevant in the modeling of the remarkable features
of equatorial flows. One such occurrence is the Equatorial Undercurrent (EUC),
a current that runs throughout the extent of the equatorial Pacific. While the
surface flow is–due to the prevailing trade-winds– preponderantly westward, a
flow reversal occurs beneath the surface leading to an eastward-flowing jet
whose core resides more of less on the thermocline and which is confined to
depths of $100-200$ m, cf. [11, 35, 36, 49]. Deeper down, at depths in excess
of about $240$ m, there is, essentially, an abyssal layer of still water.

Let us denote with $u_{w}$ the magnitude of the westward surface velocity and
with $u_{e}$ the maximal eastward velocity at the core of the EUC jet. A
simple choice for $u(r,\pi/2)$ and $u_{1}(r,\pi/2)$ that accommodates the
previously described traits of EUC is the parabolic profile

$u_{1}(r,\pi/2)=u_{e}-(u_{e}+u_{w})\left(\frac{r-R_{1}}{R_{0}-R_{1}}\right)^{2}\,\,{\rm
for}\,\,R_{1}\leq r\leq R_{0},$ (32)

and

$u(r,\pi/2)=\left\\{\begin{array}[]{ll}u_{e}-(u_{e}+u_{w})\left(\frac{r-R_{1}}{R_{0}-R_{1}}\right)^{2}&{\rm
for}\quad\overline{R}\leq r\leq R_{1}\\\ 0&{\rm for}\quad
r\leq\overline{R}\end{array}.\right.$ (33)

Then, setting

$F_{1}(r)=\sqrt{\rho_{1}(r)}\left[\Omega
r+u_{e}-(u_{e}+u_{w})\left(\frac{r-R_{1}}{R_{0}-R_{1}}\right)^{2}\right]\,\,{\rm
for}\,\,R_{1}\leq r\leq R_{0},$ (34)

and

$F(r)=\left\\{\begin{array}[]{ll}\sqrt{\rho(r)}\left[\Omega
r+u_{e}-(u_{e}+u_{w})\left(\frac{r-R_{1}}{R_{0}-R_{1}}\right)^{2}\right]&{\rm
for}\quad\overline{R}\leq r\leq R_{1},\\\ \sqrt{\rho(r)}\Omega r&{\rm
for}\quad r\leq\overline{R}\end{array}\right.$ (35)

leads to an elementary model of EUC. Here $r=R_{0}$ denotes the free surface
at the Equator, $r=R_{1}$ is the depth at which EUC attains its maximum speed
(typically this coincides with the thermocline) and

$\overline{R}:=R_{1}-(R_{0}-R_{1})\sqrt{\frac{u_{e}}{u_{e}+u_{w}}}$

represents the depth below which the fluid motion ceases, cf. Figure 2.

Figure 2: A schematic displaying the salient features of typical equatorial
current profiles in a band of about $2^{\circ}$ latitude around the Pacific
equator: a significant fluid stratification, resulting in a pycnocline (also
called thermocline) and a depth dependent current that changes from a westward
flow near the surface to a quite significant flow to the east, located below
the surface. The fluid motion dies out at depths beneath 240 m.

##  4 Analysis of the exact solutions

This section is concerned with an in-depth analysis of the exact solutions
derived previously. At first, we would like to reinforce the relevance of the
cyclostrophic balance (15), which established the connection between the
surface defining function and the pressure applied on the free surface. The
next theorem allows us to reasonably conclude that our exact solution in a
rotating frame with surface boundary conditions displays generally acceptable
and expected physical properties. More precisely, it will be demonstrated that
a growth in pressure generates a decline in height of the free surface away
from the Equator.

######  Theorem 4.1 (Monotonicity properties)

Let $\mathcal{k},\mathcal{P}_{1}$ denote the shape of the free surface and the
pressure applied on the free surface, respectively. If
$\theta_{0}\in\left(\frac{\pi}{2},\frac{\pi}{2}+\varepsilon\right)$ denotes
the polar angle of some location away from the Equator, then

$\mathcal{P}_{1}^{\prime}(\theta)<0\quad\textrm{if}\quad\mathcal{k}^{\prime}(\theta)\geq
0,$ (36)

and

$\mathcal{k}^{\prime}(\theta)<0\quad\textrm{if}\quad\mathcal{P}_{1}^{\prime}(\theta)\geq
0.$ (37)

### Proof:

To begin with, we notice that an iterative bootstrapping procedure [2] allows
us to transfer smoothing properties from $\mathcal{P}_{1}$ to $\mathcal{k}$.
By a differentiation with respect to $\theta$ in (20) we obtain

$\begin{split}P_{atm}\mathcal{P}_{1}^{\prime}(\theta)=&\frac{F_{1}^{2}\left((R_{0}+k(\theta))\sin\theta\right)}{(R_{0}+k(\theta))\sin\theta}\left(k^{\prime}(\theta)\sin\theta+(R_{0}+k(\theta))\cos\theta\right)\\\
\\\
&-\frac{F_{1}^{2}\left((R_{1}+h(\theta))\sin\theta\right)}{(R_{1}+h(\theta))\sin\theta}\left(h^{\prime}(\theta)\sin\theta+(R_{1}+h(\theta))\cos\theta\right)\\\
\\\
&-g\left(\rho_{1}(R_{0}+k(\theta))k^{\prime}(\theta)-\rho_{1}(R_{1}+h(\theta))h^{\prime}(\theta)\right)\\\
\\\
&+\frac{F_{1}^{2}\left((R_{1}+h(\theta))\sin\theta\right)h^{\prime}(\theta)}{R_{1}+h(\theta)}+F_{1}^{2}\left((R_{1}+h(\theta))\sin\theta\right)\cot\theta\\\
\\\ &-g\rho_{1}(R_{1}+h(\theta))h^{\prime}(\theta).\end{split}$ (38)

Hence, due to cancellations in the above formula, we obtain

$\begin{split}P_{atm}\mathcal{P}_{1}^{\prime}(\theta)=&\left(\frac{F_{1}^{2}\left((R_{0}+k(\theta))\sin\theta\right)}{R_{0}+k(\theta)}-g\rho_{1}(R_{0}+k(\theta))\right)k^{\prime}(\theta)\\\
\\\ &+F_{1}^{2}\left((R_{0}+k(\theta))\sin\theta\right)\cot\theta\\\ \\\
=&\left(\frac{\big{[}u_{1}(R_{0}+k(\theta),\theta)+\Omega(R_{0}+k(\theta))\sin\theta\big{]}^{2}}{R_{0}+k(\theta)}-g\right)\rho_{1}(R_{0}+k(\theta))k^{\prime}(\theta)\\\
\\\ &+F_{1}^{2}\left((R_{0}+k(\theta))\sin\theta\right)\cot\theta\\\ \\\
=&\left(\frac{\big{[}u_{1}(R_{0}+k(\theta),\theta)+\Omega(R_{0}+k(\theta))\sin\theta\big{]}^{2}}{1+\mathcal{k}(\theta)}-gR_{0}\right)\rho_{1}(R_{0}+k(\theta))\mathcal{k}^{\prime}(\theta)\\\
\\\ &+F_{1}^{2}\left((R_{0}+k(\theta))\sin\theta\right)\cot\theta.\end{split}$
(39)

Employing a similar type of argument as in Remark 3.1 and Remark 3.2 we notice
that

$\frac{\big{[}u_{1}(R_{0}+k(\theta),\theta)+\Omega(R_{0}+k(\theta))\sin\theta\big{]}^{2}}{1+\mathcal{k}(\theta)}-gR_{0}<0,$

for velocities $u_{1}$ which occur in realistic physical scenarios. The latter
inequality and relation (39) prove (36) and (37). $\square$

Another interesting and relevant upshot of the solutions proved to exist by
Theorem 3.3 refers to the regularity of the interface defining function.

######  Theorem 4.2 (Regularity of the interface)

Let us assume that the distributions of density $\rho$ and $\rho_{1}$ are
related by the equality $\rho=\rho_{1}(1+\sigma)$, where
$\tau\rightarrow\sigma(\tau)$ is a positive function satisfying
$\sigma=\mathcal{O}(10^{-3})$, cf. [38]. Then, the interface defining function
$h$ belongs to $C^{\infty}\left[\pi/2,\pi/2+\varepsilon\right]$.

### Proof:

We invoke equation (19) and differentiate with respect to $\theta$ in

$\mathcal{G}(\mathcal{h})(\theta)=0\,\,{\rm for}\,\,{\rm all}\,\,\theta,$

obtaining

$\displaystyle\Big{[}\frac{F^{2}((R_{1}+R_{1}\mathcal{h}(\theta))\sin\theta)-F_{1}^{2}((R_{1}+R_{1}\mathcal{h}(\theta))\sin\theta)}{1+\mathcal{h}(\theta)}$
(40) $\displaystyle-
gR_{1}\big{(}\rho(R_{1}+h(\theta))-\rho_{1}(R_{1}+h(\theta))\big{)}\Big{]}\mathcal{h}^{\prime}(\theta)$
$\displaystyle+\left[F^{2}((R_{1}+R_{1}\mathcal{h}(\theta))\sin\theta)-F_{1}^{2}((R_{1}+R_{1}\mathcal{h}(\theta))\sin\theta)\right]\cot\theta=0.$

We are going to prove that the bracket that multiplies
$\mathcal{h}^{\prime}(\theta)$ in the above formula is strictly positive. To
this end we notice first that

$\displaystyle\frac{F^{2}((R_{1}+R_{1}\mathcal{h}(\theta))\sin\theta)-F_{1}^{2}((R_{1}+R_{1}\mathcal{h}(\theta))\sin\theta)}{1+\mathcal{h}(\theta)}=$
(41)
$\displaystyle\rho(R_{1}+R_{1}\mathcal{h})\frac{\big{(}u(R_{1}+R_{1}\mathcal{h},\theta)+\Omega
R_{1}(1+\mathcal{h})\sin\theta\big{)}^{2}}{1+\mathcal{h}}$
$\displaystyle-\rho_{1}(R_{1}+R_{1}\mathcal{h})\frac{\big{(}u_{1}(R_{1}+R_{1}\mathcal{h},\theta)+\Omega
R_{1}(1+\mathcal{h})\sin\theta\big{)}^{2}}{1+\mathcal{h}}.$

Utilizing the relations between $\rho$ and $\rho_{1}$, as in the statement of
the Theorem, we see that the right hand side of (41) equals

$\displaystyle\frac{\rho_{1}}{1+\mathcal{h}}(u^{2}-u_{1}^{2}+2\Omega
R_{1}u(1+\mathcal{h})\sin\theta-2\Omega u_{1}R_{1}(1+\mathcal{h})\sin\theta)$
(42) $\displaystyle+\rho_{1}\sigma\frac{(u+\Omega
R_{1}(1+\mathcal{h})\sin\theta)^{2}}{1+\mathcal{h}}$
$\displaystyle<\frac{\rho_{1}u^{2}}{1+\mathcal{h}}+2\Omega
R_{1}\rho_{1}u\sin\theta+\rho_{1}\sigma\frac{(u+\Omega
R_{1}(1+\mathcal{h})\sin\theta)^{2}}{1+\mathcal{h}}$
$\displaystyle<R_{1}\left[\frac{\rho_{1}u^{2}}{R_{1}}+2\Omega\rho_{1}u\sin\theta+\rho_{1}\sigma\frac{(u+\Omega
R_{1}(1+\mathcal{h})\sin\theta)^{2}}{R_{1}}\right]<R_{1}\,kg\cdot m^{-1}\cdot
s^{-2},$

where, for the last inequality, we took into account the realistic physical
values of the involved quantities. On the other hand,

$gR_{1}(\rho-\rho_{1})=gR_{1}\rho_{1}\sigma>4.9R_{1}\,kg\cdot m^{-1}\cdot
s^{-2}.$ (43)

Owing to (42) and (43) we conclude that, indeed, the coefficient of
$\mathcal{h}^{\prime}(\theta)$ in (40) is strictly negative. We now see that,
since already $\mathcal{h}\in C^{1}\left[\pi/2,\pi/2+\varepsilon\right]$, we
can differentiate in (40) and infer that, in fact, $\mathcal{h}\in
C^{2}\left[\pi/2,\pi/2+\varepsilon\right]$. An iteration of the previous
proves the claim about the infinite differentiability of $h$. $\square$

###### Remark 4.3

We present now an explicit solution (in terms of the velocity field, the
pressure and the free surface and the interface) for an apriori given density.
More precisely, the scenario that we consider assumes that the density
$\bm{\rho}$ is as follows: throughout the bottom layer
$\\{(r,\theta,\varphi):R+d(\theta)\leq r\leq R_{1}+h(\theta)\\}$ the density
is constant and equals $\rho=\mathcal{O}(10^{3})\,kg\cdot m^{-3}$. For the
upper layer $\\{(r,\theta,\varphi):R_{1}+h(\theta)\leq r\leq
R_{0}+k(\theta)\\}$ we set

$\rho_{1}(r)=\rho-a_{1}r.$ (44)

We begin by specifying the velocity field, that is we particularize now the
functions $F$ and $F_{1}$ appearing in the formula (8) of the azimuthal
velocity $\mathbf{u}$. We will then solve a differential equation for the
interface defining function $h$, with the help of which we will find, via
(15), the formula for the free surface defining function $k$. More precisely,
we set

$F(x)=\Omega\sqrt{\rho}x\,\,{\rm and}\,\,F_{1}(x)=0\,\,{\rm for}\,\,{\rm
all}\,\,x>0.$ (45)

Theferore, the azimuthal velocity is, cf. (8)

${\bf u}(r,\theta)=\left\\{\begin{aligned} &u(r,\theta)=0\quad{\rm if}\quad
R+d(\theta)\leq r\leq R_{1}+h(\theta),\\\ &u_{1}(r,\theta)=-\Omega
r\sin\theta\quad{\rm if}\quad R_{1}+h(\theta)\leq r\leq
R_{0}+k(\theta).\end{aligned}\right.$ (46)

We will employ now (40) to derive the equation for the function $h(\theta)$
defining the interface. By means of the previous assumptions we obtain first
that

$F^{2}((R_{1}+R_{1}\mathcal{h}(\theta))\sin\theta)=\rho\Omega^{2}R_{1}^{2}(1+\mathcal{h}(\theta))^{2}\sin^{2}\theta$

and

$\rho(R_{1}+h(\theta))-\rho_{1}(R_{1}+h(\theta))=a_{1}(R_{1}+h(\theta))=a_{1}R_{1}(1+\mathcal{h}(\theta)).$

Remaining consistent with the proportion between $\rho$ and $\rho_{1}$
proposed in Theorem 4.2 we choose $a_{1}:=\frac{1}{R_{1}}\cdot 2\,kg\cdot
m^{-3}$. Availing of the previous considerations, we obtain from (40) that
$\mathcal{h}$ satisfies the equation

$\left[\rho\Omega^{2}\sin^{2}\theta-
ga_{1}\right]\mathcal{h}^{\prime}(\theta)+\rho\Omega^{2}(1+\mathcal{h}(\theta))\sin\theta\cos\theta=0.$
(47)

Using that $h(\pi/2)=0$ we solve the previous differential equation and obtain

$1+\mathcal{h}(\theta)=\sqrt{\frac{\rho\Omega^{2}-ga_{1}}{\rho\Omega^{2}\sin^{2}\theta-
ga_{1}}}.$ (48)

Taking into account that $\rho\Omega^{2}\approx 0.53\cdot 10^{-5}kg\cdot
m^{-3}\cdot s^{-2}$ and $ga_{1}\approx 0.31\cdot 10^{-5}kg\cdot m^{-3}\cdot
s^{-2}$ the equation of the interface is

$r=R_{0}+h(\theta)=R_{0}\sqrt{\frac{0.22}{0.53\sin^{2}\theta-0.31}}.$ (49)

The previous formula, forecasts a rise in the height of the interface from
$R_{0}$ at the Equator to $10^{-5}R_{0}$ at about $20$ km south and north of
Equator, cf. the graph in Figure 3.

Figure 3: The graph of the nondimensional function $1+\mathcal{h}(\theta)$.
The range of the polar angle $\theta$ is the interval
$[\frac{\pi}{2}-0.003141,\frac{\pi}{2}+0.003141]$. This corresponds in,
physical variables, to strips of about 20 km width on each side of the
Equator.

We now proceed to the explicit determination of the pressure in the upper
layer $\Omega_{top}$. It will turn out that, for the special choices (44) and
(45) we can determine an explicit formula for the pressure in the upper layer
too. To this end, we note first that, since we chose $F_{1}=0$, we obtain from
(13) and (14) the formula

$\displaystyle p_{1}(r,\theta)=$
$\displaystyle-g\int_{R_{1}+h(\theta))}^{r}\rho_{1}(\tilde{r})d\tilde{r}-g\int_{\pi/2}^{\theta}\rho_{1}(R_{1}+h(\tilde{\theta}))h^{\prime}(\tilde{\theta})d\tilde{\theta}$
(50) $\displaystyle=$ $\displaystyle-g\int_{R_{1}+h(\theta))}^{r}(\rho-
a_{1}\tilde{r})d\tilde{r}-g\int_{R_{1}}^{R_{1}+h(\theta)}(\rho-a_{1}z)dz$
$\displaystyle=$ $\displaystyle-\rho g(r-R_{1})+\rho
gh(\theta)+\frac{a_{1}g}{2}[r^{2}-(R_{1}+h(\theta))^{2}]$ $\displaystyle-\rho
gh(\theta)+\frac{a_{1}g}{2}(R_{1}+h(\theta))^{2}-\frac{a_{1}gR_{1}^{2}}{2}$
$\displaystyle=$ $\displaystyle-\rho
g(r-R_{1})+\frac{a_{1}g}{2}(r^{2}-R_{1}^{2}).$

Now, the Bernoulli relation (15) allows us to obtain the formula for the free
surface defining function $k(\theta)$. More precisely, imposing a pressure
$P_{1}(\theta)$ on the free surface, the surface defining function
$\theta\rightarrow k(\theta)$ satisfies the equation

$\frac{ga_{1}}{2}(R_{0}+k(\theta))^{2}-g\rho(R_{0}+k(\theta))+g\rho
R_{1}-\frac{ga_{1}R_{1}^{2}}{2}-P_{1}(\theta)=0\,\,{\rm for}\,\,{\rm
all}\,\,\theta.$ (51)

We treat the previous equation as a second degree equation in the unknown
$R_{0}+k(\theta)$ whose discriminant equals $g^{2}(\rho-
a_{1}R_{1})^{2}+2ga_{1}P_{1}(\theta)$. Taking into account the reasonable
physical size range for $R_{0}+k(\theta)$ we infer that

$R_{0}+k(\theta)=\frac{g\rho-\sqrt{g^{2}(\rho-
a_{1}R_{1})^{2}+2ga_{1}P_{1}(\theta)}}{ga_{1}}$ (52)

is the only possible solution to equation (51).

##  5 Conclusions

The analysis performed here provided explicit and exact steady solutions in
spherical coordinates describing flows with a free surface and an interface
which separates fluid regions of discontinuous densities. These flows manifest
an azimuthal propagation direction with no variation in the azimuthal
direction. Nevertheless, the azimuthal component of the velocity accommodates
any arbitrary variation of speed with depth below the surface: indeed the
functions $F$ and $F_{1}$ appearing in formula (8), can be chosen arbitrary.
For any such choice, $w$ and $w_{1}$ from (8) provide a solution of the Euler
equations (6). The usefulness of these solutions is rendered by their ability
to integrate observed prominent features of remarkable oceanic flows like EUC:
evidence in support of this assertion was provided in Remark 3.4.

The results in this paper represent an advancement toward the analytical
understanding of geophysical water flows. Indeed, allowing for a discontinuous
density (that gives rise to an interface) was not addresed in previous papers
[11, 12, 13, 28, 29] which handled geophysical flows with a preferred
propagation direction.

The constraints that we impose on the solutions refer to the continuity of the
pressure along the interface and the requirement that the pressure on the
surface matches some given, but otherwise, arbitrary function. The latter
leads to the Bernoulli relation on the free surface and shows that variations
in pressure give rise to changes in the shape of the free surface. While the
azimuthal velocity and the pressure function in the lower layer
$\Omega_{bottom}$ are given explicitly, the pressure in the top layer
$\Omega_{top}$ and the interface and surface defining functions, $h$ and $k$,
respectively, appear (in general) in implicit form, cf. equations (15) and
(17). However, for some specific choice of the density function and of the
azimuthal velocity, cf. (45) and (46), we were able to provide a fully
explicit solution describing an azimuthal flow: that is, also the interface
and the free surface were given explicitly, cf. (50), (52).

One aspect that needs to be considered in future works refers to the lack of
explicit formulas for the surface and interface defining functions in the
general case. This drawback was, to some degree, alleviated in Section 4 by
providing qualitative properties that reinforce the relevance of our exact
solution. More precisely, we showed that a growth in the pressure along the
free surface gives rise to a decline in height of the latter away from the
Equator. We have also proved that the interface has infinite regularity.

Future endeavours refer to the need to address some other type of solutions
which allow for more general velocity fields (e. g. having two or all three
components non-vanishing) which are relevant in the study of other realistic
ocean flows. For instance, one might want to allow a slow evolution in the
azimuthal direction; this additional dependence on the slow variable will then
yield a three-dimensional flow.

Acknowledgements: The author gratefully acknowledges the support of the
Austrian Science Fund (FWF) through research grant P 33107-N. Comments and
suggestions from anonymous referees are gratefully acknowledged. The author
would like to thank Prof. B. Basu (School of Engineering, Trinity College
Dublin) for the technical support in producing Figure 3.

Data availability: The data that support the findings of this study are
available from the corresponding author upon reasonable request.

## References

  * [1] B. Basu. On an exact solution of a nonlinear three-dimensional model in ocean flows with equatorial undercurrent and linear variation in density. _Disc. Cont. Dyn. Sys.- Ser. A_ , 2019, 39(8), pp. 4783-4796
  * [2] M. S. Berger. Nonlinearity and Functional Analysis. Academic Press, New York, 1977.
  * [3] J. P. Boyd. _Dynamics of the Equatorial Ocean_. Spinger, Berlin 2018.
  * [4] J. Chu, D. Ionescu-Kruse and Y. Yang. Exact solution and instability for geophysical waves at arbitrary latitude. _Disc. Cont. Dyn. Syst.- Ser. A_ , 39 (2019), no. 9, 5521-5541.
  * [5] A. C. Compelli, R. I. Ivanov, C. I. Martin and M. Todorov. _Surface waves over currents and uneven bottom_. Deep-Sea Research. Part II: Topical Studies in Oceanography. 160 (2019), 25-31.
  * [6] A. Constantin. An exact solution for equatorially trapped waves. _Journal of Geophysical Research: Oceans_ , vol. 117 (2012), C05029.
  * [7] A. Constantin. Some three-dimensional nonlinear equatorial flows. _J. Phys. Oceanogr._ 43 (2013) no. 1, 165-175.
  * [8] A. Constantin and R. S. Johnson. The dynamics of waves interacting with the Equatorial Undercurrent. _Geophysical and Astrophysical Fluid Dynamics_. 109 (2015), no. 4, 311–358.
  * [9] A. Constantin and R. I. Ivanov. A Hamiltonian approach to wave-current interactions in two-layer fluids. _Physics of Fluids_ 27 (2015), 086603.
  * [10] A. Constantin, R. I. Ivanov and C. I. Martin. Hamiltonian formulation for wave-current interactions in stratified rotational flows. _Arch. Ration. Mech. Anal._ 221 (2016), 1417–1447.
  * [11] A. Constantin and R. S. Johnson. An exact, steady, purely azimuthal equatorial flow with a free surface. _J. Phys. Oceanogr_. 46 (2016), no. 6, 1935-1945.
  * [12] A. Constantin and R. S. Johnson. An exact, steady, purely azimuthal flow as a model for the Antarctic Circumpolar Current. _J. Phys. Oceanogr._ 46 (2016), 3585–3594.
  * [13] A. Constantin and R. S. Johnson. A nonlinear, three-dimensional model for ocean flows, motivated by some observations of the Pacific Equatorial Undercurrent and thermocline. _Phys. Fluids_ 29 (2017), 056604.
  * [14] A. Constantin and R. S. Johnson. Large gyres as a shallow-water asymptotic solution of Euler’s equation in spherical coordinates. _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 473 (2017) no. 2200, 20170063.
  * [15] A. Constantin and R. I. Ivanov. Equatorial wave-current interactions. _Comm. Math. Phys._ 370 (2019) no. 1, 1-48.
  * [16] A. Constantin and R. S. Johnson. Ekman-type solutions for shallow-water flows on a rotating sphere: A new perspective on a classical problem. _Phys. Fluids_. 31 (2019) no. 2, 021401.
  * [17] J. Cullen and R. I. Ivanov. On the intermediate long wave propagation for internal waves in the presence of currents. _Eur. J. Mech. B Fluids_ 84 (2020), 325-333.
  * [18] J. Escher, A.-V. Matioc and B.-V. Matioc. On stratified steady periodic water waves with linear density distribution and stagnation points. _J. Diff. Equations._ 251 (2011), 2932–2949.
  * [19] J. Escher, P. Knopf, C. Lienstromberg and B.-V. Matioc. Stratified periodic water waves with singular density gradients. _Ann. Mat. Pura Appl._ 199 (2020), 1923-1959.
  * [20] A. V. Fedorov and J. N. Brown. Equatorial waves. In Encyclopedia of ocean sciences, edited by J. Steele, (Academic Press: New York, 2009), 3679–3695.
  * [21] A. Gill. _Atmosphere–Ocean Dynamics_. Academic, New York, 1982.
  * [22] B. Cushman-Roisin and J.-M. Beckers. _Introduction to Geophysical Fluid Dynamics: Physical and Numerical Aspects_. Academic Press, Waltham, Mass., 2011.
  * [23] S. Haziot and K. Marynets. Applying the stereographic projection to modeling of the flow of the antarctic circumpolar current. _Oceanography_ , 31 (2018) no. 3, 68-75.
  * [24] D. Henry. An exact solution for equatorial geophysical water waves with an underlying current. _Eur. J. Mech. B Fluids_ 38 (2013), 18-21.
  * [25] D. Henry and B.-V. Matioc. On the existence of steady periodic capillary-gravity stratified water waves. _Ann. Sc. Norm. Super. Pisa Cl. Sci._ (5) 12 (2013) no. 4, 955–974.
  * [26] D. Henry and A.-V. Matioc. Global bifurcation of capillary-gravity-stratified water waves. _Proc. Roy. Soc. Edinburgh Sect. A_ 144 (2014) no. 4, 775-786.
  * [27] D. Henry and C.I. Martin. Exact, purely azimuthal stratified equatorial flows in cylindrical coordinates. _Dyn. Partial Differ. Equ._ 15 (2018), no. 4, 337-349.
  * [28] D. Henry and C.I. Martin. Free-surface, purely azimuthal equatorial flows in spherical coordinates with stratification _J. Differential Equations_ 266 (2019) no. 10, 6788-6808.
  * [29] D. Henry and C. I. Martin. Azimuthal equatorial flows with variable density in spherical coordinate. _Arch. Ration. Mech. Anal._ 233 (2019) 497-512.
  * [30] D. Henry and C. I. Martin. Stratified equatorial flows in cylindrical coordinates. _Nonlinearity_ 33 (2020) no. 8, 3889–3904.
  * [31] H.-C. Hsu. An exact solution for nonlinear internal equatorial waves in the f-plane approximation. _J. Math. Fluid Mech._ 16 (2014), no. 3, 463-471.
  * [32] H.-C. Hsu. An exact solution for equatorial waves. _Monatsh. Math._ 176 (2015), no. 1, 143-152.
  * [33] D. Ionescu-Kruse. Exponential profiles producing genuine three-dimensional nonlinear flows relevant for equatorial ocean dynamics. _Journal of Differential Equations_ 268 (2020) no. 4, 1326-1344.
  * [34] R. I. Ivanov. Hamiltonian model for coupled surface and internal waves in the presence of currents. _Nonlinear Analysis: Real World Applications_. 34 (2017), 316-334.
  * [35] T. Izumo. The equatorial current, meridional overturning circulation, and their roles in mass and heat exchanges during the El Niño events in the tropical Pacific Ocean. _Ocean Dyn_., 55 (2005), pp. 110-123.
  * [36] G.C. Johnson, M.J. McPhaden, E. Firing. Equatorial Pacific ocean horizontal velocity, divergence, and upwelling. _J. Phys. Oceanogr_., 31 (2001), 839-849.
  * [37] R. S. Johnson. Applications of the ideas and techniques of classical fluid mechanics to some problems in physical oceanography. _Phil. Trans. R. Soc. A Math. Phys. Eng. Sci._ 376 (2018), 20170092.
  * [38] W. S. Kessler and M. J. McPhaden. Oceanic equatorial waves and the 1991-1993 El Niño. _J. Clim._ 8 (1995), 1757-1774.
  * [39] T. Lyons. Geophysical internal equatorial waves of extreme form._Discrete Contin. Dyn. Syst. Ser. A._ 39 (2019) no. 8, 4471-4486.
  * [40] J. Marshall and R. A. Plumb, Atmosphere, Ocean and Climate Dynamics: An Introductory Text, Academic, New York, 2016.
  * [41] C. I. Martin and R. Quirchmayr. Explicit and exact solutions concerning the Antarctic Circumpolar Current with variable density in spherical coordinates. _J. Math. Phys._ 60 (2019) no. 10, 101505.
  * [42] C. I. Martin and R. Quirchmayr. A steady stratified purely azimuthal flow representing the Antarctic Circumpolar Current. _Monatshefte für Mathematik_ , 192 (2020), no. 2, 401-407.
  * [43] A.-V. Matioc. An exact solution for geophysical equatorial edge waves over a sloping beach, _J. Phys. A_ 45 (2012), no. 36, 365501, 10 pp.
  * [44] A.-V. Matioc. Exact geophysical waves in stratified fluids, _Appl. Anal._ 92 (2013), no. 11, 2254-2261.
  * [45] A.-V. Matioc. On the particle motion in geophysical deep water waves traveling over uniform currents, _Quart. Appl. Math._ 72 (2014), no. 3, 455-469.
  * [46] B.-V. Matioc. Global bifurcation for water waves with capillary effects and constant vorticity. _Monatsh. Math._ 174 (2014) no. 3, 459–475.
  * [47] J. P. McCreary. Modeling equatorial ocean circulation. _Ann. Rev. Fluid Mech._ 17 (1985), 359-409.
  * [48] R. Quirchmayr. A steady, purely azimuthal flow model for the Antarctic Circumpolar Current. _Monatshefte für Mathematik_ 187 (2018) no. 3, 565-572.
  * [49] J. Sirven. The equatorial undercurrent in a two layer shallow water model. _J. Mar. Syst_., 9 (1996), 171-186.
  * [50] G. Vallis, Atmospheric and Oceanic Fluid Dynamics: Fundamentals and Large-Scale Circulation, Cambridge University Press, 2017.
  * [51] S. Walsh. Stratified steady periodic water waves. _SIAM J. Math. Anal._ 41 (2009), 1054–1105.
  * [52] S. Walsh. Steady stratified periodic gravity waves with surface tension II: global bifurcation. _Discrete Contin. Dyn. Syst._ 34 (2014) no. 8, 3287–3315.