∎

11institutetext: Jesús M. Rodríguez-de-Vera and José M. García and Gregorio
Bernabé 22institutetext: Computer Engineering Department. University of
Murcia, Murcia (Spain)
22email<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>33institutetext:
Josefa González-Carrillo 44institutetext: Hospital Virgen de la Arrixaca,
Murcia (Spain)
44email<EMAIL_ADDRESS>
# Deep learning approach to left ventricular non-compaction measurement

Jesús M. Rodríguez-de-Vera  Josefa González-Carrillo  José M. García  Gregorio
Bernabé

(Received: date / Accepted: date)

###### Abstract

Left ventricular non-compaction (LVNC) is a rare cardiomyopathy characterized
by abnormal trabeculations in the left ventricle cavity. Although traditional
computer vision approaches exist for LVNC diagnosis, deep learning-based tools
could not be found in the literature. In this paper, a first approach using
convolutional neural networks (CNNs) is presented. Four CNNs are trained to
automatically segment the compacted and trabecular areas of the left ventricle
for a population of patients diagnosed with Hypertrophic cardiomyopathy.
Inference results confirm that deep learning-based approaches can achieve
excellent results in the diagnosis and measurement of LVNC. The two best CNNs
(U-Net and Efficient U-Net B1) perform image segmentation in less than
$0.2\text{\,}\mathrm{s}$ on a CPU and in less than $0.01\text{\,}\mathrm{s}$
on a GPU. Additionally, a subjective evaluation of the output images with the
identified zones is performed by expert cardiologists, with a perfect visual
agreement for all the slices, outperforming already existing automatic tools.

###### Keywords:

Left ventricular non-compaction Medical image segmentation Cardiac image
segmentation Deep learning Convolutional neural network

##  1 Introduction

Left ventricular non-compaction (LVNC) is a recently classified form of
cardiomyopathy, characterized by abnormal trabeculations in the left ventricle
cavity Towbin et al. [2015]. It is a rare condition associated with different
types of heart defects: in Stähli et al. [2013], 12% of patients with LVNC had
other congenital cardiac diseases. In Biagini et al. [2006], 93% of patients
diagnosed with ventricular non-compaction had another cardiomyopathy.

There is still no consensus in the medical community on the way to quantify
and value LVNC Towbin et al. [2015]. Some of the already existing approaches
are the measurement of the trabecular mass based on endocardium delineation
described in Jacquier et al. [2010], or the estimation of endocardial borders
complexity using fractal methods proposed in Captur et al. [2014, 2013]. In
Bernabé et al. [2018, 2015], an automatic software tool based on medical
experience was developed to quantify the degree of non-compaction using
cardiac MRI. This tool uses traditional computer vision methods and requires
manual tuning by the cardiologists. In Bernabé et al. [2018], a self-
optimization software tool was proposed, achieving no diagnostically
significant differences in 81.34% of the images evaluated.

The irruption of deep learning (DL) techniques in the last years has led to
great improvements in medical image analysis, surpassing existing state-of-
the-art approaches Litjens et al. [2017]. Particularly, cardiac MRI
segmentation using DL has attracted considerable interest, with a growing
number of publications on the subject Chen et al. [2020]. Nevertheless,
existing papers focus mostly on finding main heart structures such as
ventricles and atria (where DL methods have equalled and even surpassed expert
analysis Bernard et al. [2018]). There are fewer works on segmenting abnormal
cardiac tissue regions which are related to some heart diseases, such as
trabeculae and LVNC.

In this work, we present a DL-based approach to the measurement and diagnosis
of LVNC. This project aims to take advantage of the significant progress in
medical image segmentation to the estimation and diagnosis of LVNC. To this
end, four different convolutional neural networks (CNNs) have been trained,
and they have been evaluated in terms of its segmentation accuracy and the
usefulness of their output to cardiologists. To the best of our knowledge,
this is the first attempt to apply DL to the problem of measuring hyper-
trabeculation levels, so we present this paper as a possible starting point
for an accurate, reliable, and fully automated diagnosis of this
cardiomyopathy.

Section 2 presents different deep learning techniques already used for image
segmentation generally, and particularly for cardiac segmentation. Section 3
describes data preprocessing, and the neural network architectures and
training methods applied to this project. The obtained performance of the
proposed method, along with its segmentation and diagnosis ability, can be
found in Section 4, where our solution is also compared with already existing
techniques. Finally, Section 5 summarizes the work, presents the main
conclusions of the paper, and proposes a future path of research.

##  2 Background

One of the most critical elements of any training process is the data fed to
the model. In particular, most of the current challenges of medical image
analysis are related to data itself Chen et al. [2020]; Hesamian et al.
[2019]; Litjens et al. [2017]. Therefore, the first element of the image
processing pipeline is usually focused on overcoming these recurrent problems.
When talking about medical image segmentation, the scarcity of annotated
images makes it difficult to create models with generalization capability, and
thus, overfitting occurs Chen et al. [2020]; Hesamian et al. [2019]. Multiple
approaches to mitigate this difficulty are found in the literature, being pre-
processing and data augmentation the most prominent ones Litjens et al.
[2017].

Given that the most usual approach to cardiac image segmentation is semantic
segmentation, the most widely used loss functions treat the problem as a
pixel-wise classification task like cross-entropy loss. For example, it is
usual to use loss functions that are intended to optimize the Dice
coefficient, like the generalized Dice loss Sudre et al. [2017] and the
Lovász-Softmax loss, presented in Berman et al. [2018]. Besides being a
differentiable surrogate function of Jaccard index, minimization of Lovász-
Softmax loss also leads to Dice coefficients optimization, achieving better
results than Dice loss in some medical image segmentation tasks Bertels et al.
[2019].

On the other hand, there also exist loss functions based on borders
differences instead of regions Kervadec et al. [2019] (boundary losses), which
make network output to be more anatomically suitable Kervadec et al. [2019].
That is important because even state-of-the-art cardiac segmentation solutions
produce anatomically impossible results. This problem is more frequent in
apical and basal slices than in mid-cavity images Bernard et al. [2018].
Moreover, the boundary losses also help to prevent the class imbalance
problem, which is particularly critical since that small tissues are usually
essential for the diagnosis.

CNNs constitute state-of-the-art for medical image segmentation (and cardiac
MRI segmentation in particular). One of the first architectures specifically
designed for medical image segmentation is the well-known U-Net Ronneberger et
al. [2015]. This network has an encoder-decoder structure, in which skip
connections between encoder and decoder are used to recover spatial context
loss and provide more precise segmentation. There are several variations of
U-Net, some of which obtain state-of-the-art results Chen et al. [2020], such
as Attention U-Net, which combines a U-shaped architecture with attention
gates Oktay et al. [2018]. Other modifications of vanilla U-Net pretend to
improve its segmentation ability or efficiency by replacing its backbone
(encoder part).

##  3 Convolutional Neural Networks for LVNC measurement

###  3.1 Network architectures

Given that trabecular measurement using neural networks is an issue that has
never been addressed, it is unclear which architecture would achieve the best
performance. For that reason, this work proposes the usage of four different
network architectures in the experiments, whose results are compared to
identify which one performs better for this problem. These four architectures
are as follows:

  * •

The vanilla U-Net Ronneberger et al. [2015], as it is widely used for cardiac
image segmentation Chen et al. [2020] and it is the basis of many state-of-
the-art solutions for medical segmentation.

  * •

Attention U-Net Oktay et al. [2018]. It includes attention gates into the
U-Net architecture. Attention gates help the model to focus on the most
important areas of the image by reducing the contribution of non-relevant
regions.

  * •

Efficient U-Net B1, which uses the architecture EfficientNet-B1 Tan and Le
[2019] as a substitute for the encoder part of classical U-Net. This network
is included as a lightweight architecture alternative to the other models.

  * •

Pyramid Attention Network (PAN) Li et al. [2018]. It is a lightweight encoder-
decoder network, which uses attention mechanisms (similar to the attention
gates but more lightweight) and presents an alternative to spatial pyramid
pooling. We included this network because this kind of blocks (attention and
pyramid pooling) have been used in cardiac image segmentation, achieving
state-of-the-art results in some cases Chen et al. [2020].

In order to estimate and compare models goodness, three kind of evaluations
metrics have been utilized: a) Inference time when running both in GPU and
CPU, b) Segmentation capabilities (based on Dice coefficient), and c)
Diagnosis ability based on network output.

It is important to notice that none of the used architectures takes into
account temporal or spatial context.

###  3.2 Non-Compaction Measurement and LVNC Diagnosis

The principal purpose of this work is to diagnosis LVNC in a fully automatic
way, based on available data, i.e. heart images (more specifically MRI
slices). In order to measure LVNC levels, three different regions of the left
ventricle must be distinguished (as shown in Figure 1): external layer (EL),
internal cavity (IC) and trabeculae (T).

Figure 1: Segmentation sample of left ventricle. Green: external layer.
Yellow: trabeculae or non-compacted tissue. Blue: internal cavity

The area of these parts is then used to compute the percentage quantification
of the trabecular zone with respect to the compacted area with the following
expression Bernabé et al. [2018]:

$PTA=100\cdot\frac{TA}{TA+ELA}\quad[\%]$ (1)

where:

$TA=\text{trabeculae area}\quad\text{and}\quad ELA$ ${}={}$ | external layer area
---|---

Although IC information does not appear explicitly in the above expression, we
included it as a separate class to distinguish it from the background.

Once the value of expression (1) is known for an image, the slice is
classified using a threshold: values greater than or equal to 27.4% are
considered as a LVNC indicative Bernabé et al. [2017].

###  3.3 Data

The original data set was composed of 2857 heart short-axis MRI slices of 292
patients, which were stored in DICOM files. All cases had previously been
diagnosed with hypertrophic cardiomyopathy (HCM) which is another
cardiomyopathy in which myocardium becomes abnormally thick.

Magnetic resonance studies were performed in two hospitals with different
scanners:

  * •

Magnetic Resonance of the Southeast. Affidea (Spain). Imaging was done on a
1.5T scanner (SIGNA HDxt 1,5: General Electric Systems, USA). Image dimension
is $224\times 224$ for every subject and pixel spacing of $1.75\times
1.75\times 0.8mm$ with $2mm$ between slices.

  * •

University Hospital Virgen de la Arrixaca in Murcia (Spain). Imaging was done
on a 1.5T magnet (Achieva CV, Philips Medical Systems, Netherlands). Image
dimension is $256\times 256$ for every subject and pixel spacing of $1.5\times
1.5\times 0.8mm$ with $2mm$ between slices.

Studies in both hospitals were done without contrast agent, in synchronization
with the ECG and patients were holding their breath (at the end of an exhale).
Moreover, slices generated with the first scanner were resampled to a
resolution of $256\times 256$, as it corresponds to the smallest pixel
spacing.

Given the importance of the data set in neural networks training, image slices
were successively filtered in a process which we describe in the next lines.

Firstly, we discarded slices where left ventricle was not visible, which
resulted in the reduction of the data set to 2438 slices.

Then, we proceeded to the annotation step. In this research, no traditional
and time expensive fully manual annotation was performed. Instead of that,
images were labeled using the tool QLVTHC developed in Bernabé et al. [2018]
to segment the whole data set, and cardiologists just manual tuned its
parameters for each slice. By doing this, we pretend to overcome a typical
issue of this field: this kind of projects usually requires a human expert to
segment the whole data set manually. The segmentation produced by this tool
has a resolution of $512\times 512$.

After that, each segmentation produced by QLVTHC, which contained the three
regions showed in Figure 1, was evaluated by two cardiologists using the 1 to
5 scale proposed in Gibson et al. [2004] that is also used in Bernabé et al.
[2018]. This scale allows to measure the quality of the segmentation regarding
diagnosis. Thus, a value higher than or equal to $4.0$ indicates that the
segmentation does not represent diagnostically significant differences, while
values lower than or equal to $3.5$ indicate that it does. Finally, only
images with a score higher than or equal to 4.5 were kept.

Once the selection based on cardiologists evaluation was performed, the
dataset was reduced to 2362 images.

The last stage of the filtering process aims to reduce a potential problem
derived from resampling: resizing ground-truth segmentation (original size of
$512\times 512$) could lead into least sharp class maps. In order to quantify
the deformation degree of annotated segmentation, two criteria were applied:

  * •

For each slice, we computed the percentage of the external layer, the
percentage of the trabeculae and the expression (1) before and after the
resampling. Images for which the relative error of the resampled segmentation
was above $5\%$ for any of these percentages were discarded.

  * •

As a measure of topological correctness, we computed the number of connected
components of the trabecular area before and after the resampling. Only slices
for which this value did not differ were kept.

These two criteria were met by 2100 MRI slices corresponding to 277 patients,
which conformed the final data set used for training. A total of 977 slices
corresponded to patients suffering from LVNC.

In addition to the above explained image filtering, another one was performed
after the first training. More detailed explanation about it can be found in
subsection 4.2.

###  3.4 Training

We considered whole-heart models, i.e. separated models for apical, mid-cavity
and basal slices were not trained. The resampled MRI slices are normalized to
have zero mean and unit variance before they are fed to each network.
Moreover, data augmentation was performed for the training set by applying
random rotations of 90, 180 and 270 degrees with a $0.25$ probability.

The loss function which was used in all experiments consists of a linear
combination of two components:

$\mathcal{L}=2\cdot\mathcal{L}_{L}+\mathcal{L}_{B}$ (2)

where

$\mathcal{L}_{L}$ | is the Lovász-Softmax loss defined in Berman et al. [2018].
---|---
$\mathcal{L}_{B}$ | is a boundary loss inspired in Kervadec et al. [2019]. It is computed only for trabeculae.

With this expression of $\mathcal{L}$, we pretend to obtain a soft loss
function to speed up and facilitate learning while focussing on the
minimization of Sørensen–Dice coefficient (by using $\mathcal{L}_{L}$). With
$\mathcal{L}_{B}$ we intend to help addressing the class imbalance problem and
to make the output of the networks more anatomically plausible.

Because of the presence of several layers using ReLU activation function in
all of the proposed architectures, we use Kaiming initialization to overcome
gradient vanishing He et al. [2015]. For minimizing the mentioned loss
function, Radam optimizer was used, with an initial learning rate of $0.001$
and a weight decay of $0.0005$. Finally, a batch size of 24 was used for every
training.

##  4 Results and Discussion

The evaluation platform is equipped with a double socket CPU Intel Xeon
E5-2603 v3 1.60GHz (12 cores and 64 GB RAM), and two different NVIDIA GPUs:
the GeForce RTX 2080 Ti (11GB RAM, 1.350 GHz) for the training, and the
GeForce GTX 1080 (8 GB RAM, 1,607 GHz) for the inference.

###  4.1 Network Training

As has been mentioned in subsection 3.1, four different architectures were
evaluated for the problem of LVNC diagnosis. As we can see in Table 1, there
are great differences concerning the number of parameters of each
architecture. It is noteworthy that Efficient U-Net B1 has less than half
parameters than the second smallest network.

| Num. Params | Training time (min)
---|---|---
U-Net | $31,383,876$ | $42$
Attention U-Net | $34,877,616$ | $45$
Efficient U-Net B1 | $8,344,804$ | $26$
PAN | $21,470,411$ | $28$
Table 1: Number of parameters of each model and average training time in
minutes using a GPU Nvidia GeForce RTX 2080 Ti.

For training, a 5-fold cross-validation process was used, so the data set was
split into five equal-sized folds or subsets in a way that all of them had the
same proportion of slices of patients with LVNC (stratified cross-validation).
The same folds were used for every architecture. Thus, for each neural
network, we obtain five models to evaluate the architecture: one fold is used
as a test set, and the model is trained using the other four-folds (which are
splited in train-validation with 80%-20%), and this process is repeated for
each fold.

Once the model had been trained, it was evaluated over the test set (the left-
out fold). Thus, the results discussed below correspond to the whole data set,
but the output for each image was produced by a model that had never seen it
before.

The training strategy was the same for every fold and architecture: a maximum
of 25 epochs is established, but an early stopping can be triggered if loss
function $\mathcal{L}$ does not improve during 5 epochs. None of the training
processes reached the twenty-fifth epoch, obtaining a maximum of 22 epochs for
vanilla U-Net. Due to disparities in the number of parameters, we also found
variations in training time for different architectures, as shown in Table 1,
which contains (in the second column) the average training time of the five
training processes for each network.

###  4.2 First Results and Data Refinement

After the first 5-fold training of every model, outputs were analyzed. It was
found that generally, segmentations produced by the trained neural networks
were more anatomically plausible than ground truth generated by QLVTHC.
Moreover, for specific slices, the output of each of the proposed models was
significantly more accurate than reference segmentations (example in Figure
2). While it is true that these deficiencies in ground truth segmentations are
not relevant from the diagnosis point of view, the presence poor quality data
could potentially add some bias and prevent the networks from learning and
generalizing correctly. Therefore, it was hypothesized that discarding images
for which the output of trained models was better than ground truth could lead
to more accurate models. A total of 80 images were identified and subtracted
from the data set, and a new 5-fold training was performed for each
architecture.

Next, we evaluated the impact of removing these slices. Both models were
trained with the full data set (composed of 2100 MRI slices), and the reduced
data set (2020 images). The evaluation was carried out only over the reduced
one to obtain a fair comparison. Figure 3 shows the values of the Dice
coefficient for each of the regions considered in the segmentation task. The
obtained Dice score is higher when training without the removed slices for
every architecture and the three regions. For this reason, the rest of the
results and discussion in this section deals with models trained with the
reduced data set.

(a) MRI slice

(b) Ground truth

(c) U-Net output

Figure 2: Model outputs for a sample slice Figure 3: Dice coefficient of EL,
IC and T regions of left ventricle for every architecture, trained with both
the full data set and the reduced data set.

###  4.3 Inference Time and Segmentation Results

Although most of the papers about neural networks mainly focus on the time
taken by the training process, we find the measurement inference time equally
relevant. The reason for this is that the main purpose of this research is to
provide cardiologists with an automated tool that they could use with the
available equipment and that helps them by saving time (also when comparing
with existing tools). Therefore, the time required for inference was measured
when running over CPU and GPU, and it can be found in Table 2: the first
column contain results for the Intel CPU, and the second one for the NVIDIA
GPU GeForce GTX 1080.

Both time metrics were obtained from 100 runs with a previous warm-up process
of 5 runs. Each one of the architectures performs inference in less than one
second when using CPU, but we can highlight both PAN and Efficient U-Net B1
whose inference time is below $0.2\text{\,}\mathrm{s}$. Focusing on the GPU
results, the benefits are even more evident: thanks to batch processing, a
typical MRI study of 15 slices could be completed in less than
$0.05\text{\,}\mathrm{s}$ in any of the networks.

In order to better understand the impact of results in Table 2 in time savings
for cardiologists, a comparison with existing techniques is performed. Three
methods are compared with the proposed models for the already mentioned
problem of analyzing 15 slices: fully manual segmentation, QLVTHC semi-
automatic tool, and the fully automated tool proposed in Bernabé et al.
[2018], SOST. Times for fully manual and semi-automatic approaches (25 minutes
and 10 seconds per slice, respectively) have been reported by cardiologists.
The execution time of SOST has been collected from Bernabé et al. [2018], and
we picked up the lowest times declared in the paper to achieve a fairer
comparison, given that the same machine could not be used. Execution times of
both traditional tools correspond to CPU execution, as GPU version of them
does not exist.

As we can see in Figure 4, GPU execution of our models clearly outperforms the
three methods, with speedups of up to six orders of magnitude with respect to
traditional manual segmentation, and up to $612\text{\,}\mathrm{w}$hen
comparing with SOST. In regards with CPU executions, our models also provide
great time saving with respect to manual segmentation (up to
$12,068\text{\,}\mathrm{t}$imes faster), and the two more lighweight models
(PAN and Efficient U-Net B1) are faster than the existing automatic solution
SOST. Therefore, these results reinforce the idea that applying accurate and
reliable fully automatic tool would lead to a significant time saving for
specialists.

| Time CPU (ms) | Time GPU (ms)
---|---|---
U-Net | $804.6\pm 7.2$ | $7.0\pm 0.1$
Attention U-Net | $779.8\pm 28.7$ | $12.7\pm 0.1$
Efficient U-Net B1 | $198.4\pm 15.4$ | $42.9\pm 1.3$
PAN | $124.3\pm 16.8$ | $18.6\pm 0.1$
Table 2: Mean and standard deviation of inference time for a single slice on
CPU, and same metrics for batches of 24 slices on GPU. Figure 4: Estimated
speedups for the problem of quantifying PTA of 15 slices with regard to manual
annotation and both manually configured tool (QLVTHC) and fully automated tool
(SOST Bernabé et al. [2018]). Dashed line indicates unit speedup.

In Table 3, metrics about the quality of segmentations of every model are
displayed. The Dice coefficients of the three U-Net based architectures are
very similar, and it is only worthy of note the differences in the error of
computed PTA. On the other hand, the difference between the three U-Net based
models and PAN is more significant. According to Dice scores, a non-compacted
area (trabeculae) seems to be the most difficult to identify by each one of
the models. While these differences could be due to class imbalance, given
that loss function was chosen to prevent this issue, we are inclined to
believe that variations in these values result from the difference in the
intrinsic difficulty of recognizing each region.

| PTA Error (%) | Dice EL | Dice IC | Dice T
---|---|---|---|---
U-Net | $4.20\pm 4.69$ | $0.89\pm 0.11$ | $0.96\pm 0.13$ | $0.84\pm 0.16$
Attention U-Net | $4.45\pm 5.54$ | $0.89\pm 0.11$ | $0.95\pm 0.10$ | $0.84\pm 0.15$
Efficient U-Net B1 | $4.31\pm 5.17$ | $0.89\pm 0.10$ | $0.94\pm 0.09$ | $0.83\pm 0.15$
PAN | $4.65\pm 5.63$ | $0.86\pm 0.11$ | $0.91\pm 0.10$ | $0.76\pm 0.15$
Table 3: Mean metrics ($\pm$ standard deviation) of every model. Error
committed when computing PTA based on each network output and Dice score for
external layer (EL), internal cavity (IC) and trabeculae (T).

###  4.4 Diagnosis

Once a slice is segmented, the presence of LVNC is predicted based on the rule
of 27.4% threshold for PTA. Table 4 contains some metrics which are usually
used to compare diagnosis models for every architecture taken into account.
Recall and specificity are included to provide a better understanding of each
model’s performance regarding positives and negatives separately. On the other
hand, accuracy is calculated because it is easy to understand, and it is
widely used in literature; however, Matthew coefficient is used as the global
comparison metric of the four architectures given that it is more informative
than the F1 score and accuracy Chicco and Jurman [2020]. U-Net outperforms the
three others concerning accuracy and Matthews coefficient, and its specificity
and recall are slightly lower than the best values. Attention U-Net achieves
the best specificity value of the four models, but it also obtains a recall
significantly lower than the other three architectures. For its part,
Efficient U-Net B1 attains the highest recall and presents the second-best
accuracy and Matthews coefficient right behind U-Net. Finally, PAN got the
worst results in three of the four metrics, but its recall is very similar to
the best one.

| Accuracy |  | Matthews
---
Coefficient
Recall | Specificity
U-Net | $0.865$ | $0.727$ | $0.836$ | $0.889$
Attention U-Net | $0.850$ | $0.698$ | $0.797$ | $0.894$
Efficient U-Net B1 | $0.852$ | $0.703$ | $0.841$ | $0.862$
PAN | $0.847$ | $0.691$ | $0.839$ | $0.853$
Table 4: Accuracy, Matthews coefficient, recall and specificity of diagnosis
models based on each neural network using the rule of 27.4% PTA threshold.

Because of the results showed in Table 3 and Table 4, U-Net presents the best
performance because it achieves the best average Dice coefficient for each
region and the best Matthews coefficient as a diagnosis model. It is also
noteworthy that, despite Efficient U-Net B1 is substantially less complex than
the others in terms of the number of parameters, its performance is close to
U-Net and Attention U-Net. Moreover, Efficient U-Net B1 inference on CPU is
more than four times faster than U-Net, and it is even faster than traditional
automated methods, so we consider that this architecture should be taken into
account as an efficient alternative to U-Net.

Therefore, outputs of both U-Net and Efficient U-Net B1 were evaluated by two
cardiologists in order to determine the validity of the proposed models
clinically. The evaluation was performed using the scale previously commented
in subsection 3.3 Gibson et al. [2004]. It is important to notice that despite
differences in quantitative metrics visible in Table 3 and Table 4, both
architectures obtained the same punctuation for each one of the slices. This
fact indicates that our selection of the best two models that should be
evaluated by specialists is coherent in the sense that their outputs are
diagnostically equivalent in spite of being so different regarding the number
of parameters.

None of the outputs produced by either U-Net or Efficient U-Net B1 presented
diagnostically significant issues. More specifically, $88.84\%$ of the cases
were evaluated with the highest score ($5.0$), meaning that cardiologists
found the generated segmentation to be completely accurate. The rest of the
images obtained an score of $4.5$, so no noticeable differences were found
either.

If we take into account results obtained by SOST (reported in Bernabé et al.
[2018]), our two selected models outperform this previously existing automatic
tool. Indeed, our full-heart models achieve better cardiologists evaluation
even than apical, basal, and mid-cavity specific models of SOST. In
particular, the percentage of images whose output obtained a score equal or
higher than 4.5 was 78.22% for the best SOST model (only valid for mid-cavity
slices) and 63.16% for the full-heart SOST model; in comparison with the 100%
of the two networks selected in this paper. It is noteworthy that images
utilized in the evaluation of SOST in Bernabé et al. [2018] constitute a
subset of the data set used during this research.

##  5 Conclusions

In this paper, we have proposed four different convolutional neural networks
as a first approach to LVNC measurement using deep learning techniques. The
three U-Net alternatives performed better than PAN architecture for
segmentation and diagnosis metrics. We could remark two architectures:
firstly, vanilla U-Net, which obtained the best results in almost every
metric. On the other hand, Efficient U-Net B1, which presented similar metrics
to U-Net, and its number of parameters, training time and inference time on
CPU are significantly lower than other architectures.

U-Net and Efficient U-Net B1 were evaluated by cardiologists taking into
account the goodness of their outputs for diagnosis, i.e. a subjective
evaluation of the accuracy of the generated segmentations. Each one of the
slices was found to be perfectly accurate for diagnosis, and 88.84% of them
was evaluated as a perfect match.

Moreover, using the two finally selected models, segmentations can be obtained
in less than $0.2\text{\,}\mathrm{s}$ when using a CPU and in less than
$0.01\text{\,}\mathrm{s}$ on a GPU. These times corroborate that this kind of
tools could lead to considerable time savings of the specialists, given that
fully manual segmentation takes them about 25 minutes per slice.

Results achieved confirm the main idea that motivated this research: deep
learning based approaches can achieve excellent results in the diagnosis and
measurement of LVNC. In fact, proposed vanilla U-Net and Efficient U-Net B1
solutions achieve improvements in both execution time and quality of outputs
(based on cardiologists evaluation) when compared with already existing fully
automatic method.

In spite of these promising results, this paper presents only a first approach
to this problem, and critical challenges remain. Firstly, the trabecular area
presents the lowest Dice coefficient in every considered case. Given that the
problem we are trying to solve is the quantification of trabeculae, exporing
some options to improve this metric should be a priority. Moreover, it could
be interesting to use techniques like weakly supervised learning to work with
poor quality data, like the ground truth segmentations generated by already
existing tools.

###### Acknowledgements.

This work was supported by the Spanish MCIU and AEI, as well as European
Commission FEDER funds, under grant RTI2018-098156-B-C53.

## References

  * Berman et al. [2018] Berman M, Triki AR, Blaschko MB (2018) The lovasz-softmax loss: A tractable surrogate for the optimization of the intersection-over-union measure in neural networks. 2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition DOI 10.1109/cvpr.2018.00464
  * Bernabé et al. [2018] Bernabé G, Casanova JD, Cuenca J, González-Carrillo J (2018) A self-optimized software tool for quantifying the degree of left ventricle hyper-trabeculation. The Journal of Supercomputing 75(3):1625–1640, DOI 10.1007/s11227-018-2722-x
  * Bernabé et al. [2015] Bernabé G, Cuenca J, de Teruel PEL, et al. (2015) A software tool for the automatic quantification of the left ventricle myocardium hyper-trabeculation degree. Procedia Computer Science 51:610 – 619, DOI https://doi.org/10.1016/j.procs.2015.05.329, international Conference On Computational Science, ICCS 2015
  * Bernabé et al. [2017] Bernabé G, González-Carrillo J, Cuenca Muñoz J, et al. (2017) Performance of a new software tool for automatic quantification of left ventricular trabeculations. Revista Española de Cardiología (English Edition) 70(5):405 – 407, DOI 10.1016/j.rec.2016.07.006
  * Bernard et al. [2018] Bernard O, Lalande A, Zotti C, et al. (2018) Deep Learning Techniques for Automatic MRI Cardiac Multi-structures Segmentation and Diagnosis: Is the Problem Solved? IEEE Transactions on Medical Imaging 37(11):2514–2525, DOI 10.1109/TMI.2018.2837502
  * Bertels et al. [2019] Bertels J, Eelbode T, Berman M, et al. (2019) Optimizing the dice score and jaccard index for medical image segmentation: Theory and practice. In: Medical Image Computing and Computer Assisted Intervention – MICCAI 2019, Springer International Publishing, Cham, pp 92–100, DOI 10.1007/978-3-030-32245-8˙11
  * Biagini et al. [2006] Biagini E, Ragni L, Ferlito M, et al. (2006) Different types of cardiomyopathy associated with isolated ventricular noncompaction. The American journal of cardiology 98:821–4, DOI 10.1016/j.amjcard.2006.04.021
  * Captur et al. [2013] Captur G, Muthurangu V, Cook C, et al. (2013) Quantification of left ventricular trabeculae using fractal analysis. Journal of Cardiovascular Magnetic Resonance 15, DOI 10.1186/1532-429X-15-36
  * Captur et al. [2014] Captur G, Lopes LR, Patel V, Li C, et al. (2014) Abnormal cardiac formation in hypertrophic cardiomyopathy: fractal analysis of trabeculae and preclinical gene expression. Circulation: Cardiovascular Genetics 7(3):241–248, DOI 10.1161/CIRCGENETICS.113.000362
  * Chen et al. [2020] Chen C, Qin C, Qiu H, et al. (2020) Deep learning for cardiac image segmentation: A review. Frontiers in Cardiovascular Medicine 7, DOI 10.3389/fcvm.2020.00025
  * Chicco and Jurman [2020] Chicco D, Jurman G (2020) The advantages of the matthews correlation coefficient (mcc) over f1 score and accuracy in binary classification evaluation. BMC Genomics 21, DOI 10.1186/s12864-019-6413-7
  * Gibson et al. [2004] Gibson D, Spann M, Woolley S (2004) A wavelet-based region of interest encoder for the compression of angiogram video sequences. IEEE Transactions on Information Technology in Biomedicine 8(2):103–113, DOI 10.1109/titb.2004.826722
  * He et al. [2015] He K, Zhang X, Ren S, Sun J (2015) Delving deep into rectifiers: Surpassing human-level performance on imagenet classification. 2015 IEEE International Conference on Computer Vision (ICCV) DOI 10.1109/iccv.2015.123
  * Hesamian et al. [2019] Hesamian MH, Jia W, He X, Kennedy P (2019) Deep learning techniques for medical image segmentation: Achievements and challenges. Journal of Digital Imaging 32(4):582–596, DOI 10.1007/s10278-019-00227-x
  * Jacquier et al. [2010] Jacquier A, Thuny F, Jop B, et al. (2010) Measurement of trabeculated left ventricular mass using cardiac magnetic resonance imaging in the diagnosis of left ventricular non-compaction. European heart journal 31(9):1098–1104, DOI 10.1093/eurheartj/ehp595
  * Kervadec et al. [2019] Kervadec H, Bouchtiba J, Desrosiers C, et al. (2019) Boundary loss for highly unbalanced segmentation. In: Proceedings of The 2nd International Conference on Medical Imaging with Deep Learning, PMLR, London, United Kingdom, Proceedings of Machine Learning Research, vol 102, pp 285–296, URL http://proceedings.mlr.press/v102/kervadec19a.html
  * Li et al. [2018] Li H, Xiong P, An J, Wang L (2018) Pyramid attention network for semantic segmentation. 1805.10180
  * Litjens et al. [2017] Litjens G, Kooi T, Bejnordi BE, et al. (2017) A survey on deep learning in medical image analysis. Medical Image Analysis 42:60–88, DOI 10.1016/j.media.2017.07.005
  * Oktay et al. [2018] Oktay O, Schlemper J, Folgoc LL, et al. (2018) Attention u-net: Learning where to look for the pancreas. 1804.03999
  * Ronneberger et al. [2015] Ronneberger O, Fischer P, Brox T (2015) U-net: Convolutional networks for biomedical image segmentation. Lecture Notes in Computer Science Medical Image Computing and Computer-Assisted Intervention – MICCAI 2015 p 234–241, DOI 10.1007/978-3-319-24574-4˙28
  * Stähli et al. [2013] Stähli BE, Gebhard C, Biaggi P, et al. (2013) Left ventricular non‐compaction: Prevalence in congenital heart disease. International Journal of Cardiology 167(6):2477 – 2481, DOI 10.1016/j.ijcard.2012.05.095
  * Sudre et al. [2017] Sudre CH, Li W, Vercauteren T, et al. (2017) Generalised dice overlap as a deep learning loss function for highly unbalanced segmentations. Lecture Notes in Computer Science p 240–248, DOI 10.1007/978-3-319-67558-9˙28
  * Tan and Le [2019] Tan M, Le QV (2019) Efficientnet: Rethinking model scaling for convolutional neural networks. 1905.11946
  * Towbin et al. [2015] Towbin JA, Lorts A, Jefferies JL (2015) Left ventricular non-compaction cardiomyopathy. The Lancet 386:813–825, DOI 10.1016/S0140-6736(14)61282-4


# Geophysical internal equatorial waves of extreme form

<EMAIL_ADDRESS>
###### Abstract.

The existence of internal geophysical waves of extreme form is confirmed and
an explicit solution presented. The flow is confined to a layer lying above an
eastward current while the mean horizontal flow of the solutions is westward,
thus incorporating flow reversal in the fluid.

###### Key words and phrases:

Equatorial flows, geophysical internal waves, extreme waves.

###### 1991 Mathematics Subject Classification:

Primary: 76B55, 76B15; Secondary: 35Q31.

The paper is for the special theme: Mathematical Aspects of Physical
Oceanography, organized by Adrian Constantin.

Tony Lyons

Department of Computing & Mathematics

Waterford Institute of Technology

Waterford, Ireland

##  1\. Introduction

Hydrodynamic waves of extreme form have played a prominent role in the
mathematical study of water waves since the 19th century, when Stokes [39]
first predicted the existence of the wave of greatest height. Extreme waves
are inherently nonlinear, for instance extreme Stokes waves are solutions of
full Euler equations for inviscid, incompressible fluids, and their existence
(cf. [1]) relies on the nonlinear character of these governing equations.
Moreover, extreme waves incorporate numerous mathematical intricacies [3, 40],
especially in regions close to the wave crest, which presents itself in the
form of a cusp. Nevertheless, numerous physical features of such flows have
been analysed mathematically, for instance the particle trajectories and
pressure in extreme Stokes waves in deep-water and water of finite depth may
be found in [9, 28, 29, 30, 31]. In this work we will investigate internal
geophysical flows of extreme form near the equator, and it will be found that
the presence of a cusp at the wave crest of these flows once again introduces
several mathematical difficulties. In particular, the necessary criteria to
apply the implicit function theorem used to prove the existence of smooth
geophysical flows are absent in the extreme case, while the cross-section of
the wave about the equator also requires a more intricate analysis, both of
which are addressed in this work.

The mathematical analysis of geophysical waves in the equatorial region has
recently witnessed an explosion of research activity. A significant motivation
for this interest has been the fact that such waves have an exact description
in a variety of settings in the form of equatorially trapped Gerstner like
solutions, see [23] for a review of three-dimensional Gerstner like solutions
in the equatorial region. Equatorially trapped surface waves have been
investigated in [7, 13, 14] while surface waves in the presence of underlying
currents have been treated in [22, 25]. Internal geophysical flows are
investigated in [10, 4, 33, 32, 36, 27], while instabilities in such flows are
addressed in [12, 26, 18, 24].

The waves we shall consider are inherently nonlinear and while this
nonlinearity may preclude an exact analysis of their behaviour in many
instances, this is not strictly the case, with the Gerstner solution [19]
being one of the few known exact solutions of the governing equations for
incompressible, inviscid fluids in deep-water. The Gerstner wave was
subsequently rediscovered by Rankine [34] and of particular interest it was
extended to heterogeneous fluids by Dubreil-Jacotin [17]. A more modern
perspective of Gerstner’s wave is found in [5, 6, 21] where the evolution of
the fluid domain under the action of this solution is analysed using a
combination of analytical and topological methods.

Of particular significance to the current work, Gerstner’s solution may also
be extended to hydrodynamic systems described in rotating coordinate systems.
Gerstner-like solutions were initially used in [8, 10] to model three-
dimensional equatorially trapped geophysical flows and in [11] a Gerstner-like
solution was constructed as an exact solution for internal geophysical water-
waves near the equator in the $\beta$-plane formulation. This explicit
solution is constructed in terms of Lagrangian variables (see [2] for a
comprehensive review of the Lagrangian approach to fluid dynamics), and it is
found that these waves may propagate in an eastward direction and are also
found to decay significantly in the meridional direction away from the
equator. Moreover, this Gerstner-like solution exists in a layer of finite
depth, lying directly above a uniform current layer, which in turn lies above
a transitional layer where the uniform current decays to zero when the
motionless layer begins. The uniform current layer and the transitional layer
are both subject to eastward currents, while it was shown in [11] that the
Gersner-like solution induces a mean current towards the west, and so the
model accounts for flow reversal, as observed in the Pacific equatorial
undercurrent.

Notably, the results presented in [8, 10, 11] concern nonlinear flows which
are smooth throughout. In the current work we aim to extend the results
presented in [11] to incorporate the case of nonlinear extreme waves in the
Gerstner-like solution, which feature cusps at the wave-troughs along the
equator in the thermocline, the interface between the Gerstner-like solution
and the uniform current layer. The presence of these cusps leads to a number
of mathematical complications, in particular it means the criteria necessary
to confirm the existence of the extreme wave profile in the thermocline by
means of the implicit function theorem are not satisfied along the equator.
This issue is overcome in this work thus confirming the existence of extreme
waves in the thermocline. Based on this, we proceed to outline certain
features the meridional cross-section of the thermocline must display.
Finally, we also confirm that current flow reversal across the thermocline is
preserved, even when it is disturbed to the form of an extreme Gerstner-like
solution.

##  2\. Preliminaries

###  2.1. The governing equations in a rotating frame

The rotation of the Earth about the polar axis induces a fictitious force,
known as the Coriolis force, on all bodies away from the equator. In the
following, the Earth is approximated to be an ideal sphere with radius
$R=6731\,\mathrm{km}$ and rotating with angular velocity $\Omega=7.29\times
10^{-5}\,\rm{rad\,s^{-1}}$ about the polar axis. The acceleration due to
gravity is assumed uniform over Earth’s surface with numerical value
$g=9.8\,\rm{m\,s^{-2}}$.

Figure 1. The rotating $\left(x,y,z\right)$-coordinate system fixed to the
surface of the Earth. The $x$-axis points due-east, the $y$-axis points due-
north and the $z$-axis points vertically upwards from the surface.

In Figure 1 we sketch the surface of the Earth with a Cartesian coordinate
system, $(x,y,z)$ constructed at some point about the equator. The $x-$axis of
this coordinate system points eastwards parallel to the equator, while the
$y$-axis points due-north towards the pole. The $z$-axis is the line
perpendicular to the surface of the Earth at this point. The latitude of the
equator is set to be $0^{\circ}$, while the lines of latitude approximately
$\pm 5^{\circ}$ either side of the equator denote the boundary of the
equatorial region.

In the following, we shall consider a fluid in motion relative to the rotating
$(x,y,z)$-coordinate axes with associated orthonormal basis vectors
$\left\\{\mathbf{e}_{x},\mathbf{e}_{y},\mathbf{e}_{z}\right\\}$. This rotation
is defined in relation to a reference frame which is stationary with respect
to the centre of the Earth and whose orthonormal basis vectors we denote by
$\left\\{\mathbf{i},\mathbf{j},\mathbf{k}\right\\}$. The coordinate
transformation from the inertial to the rotating frame of reference is given
by

(1) $\begin{cases}\mathbf{e}_{x}=-\sin\phi\mathbf{i}+\cos\phi\mathbf{j}\\\
\mathbf{e}_{y}=-\sin\theta\cos\phi\mathbf{i}-\sin\theta\sin\phi\mathbf{j}+\cos\theta\mathbf{k}\\\
\mathbf{e}_{z}=\cos\theta\cos\phi\mathbf{i}+\cos\theta\sin\phi\mathbf{j}+\sin\theta\mathbf{k}\end{cases}$

with $\phi\in[0,2\pi]$ and
$\theta\in\left[-\frac{\pi}{2},\frac{\pi}{2}\right]$ respectively denoting the
zonal and meridional coordinates.

The inertial coordinate system may be chosen such that the rotation vector of
the Earth is given by $\mathbf{\Omega}=\Omega\mathbf{k}$, where $\mathbf{k}$
is the unit vector projecting from the north pole along the polar axis.
Inverting the transformation given by equation (1), the rotation vector
$\mathbf{\Omega}$ is given by

(2)
$\mathbf{\Omega}=\Omega\cos\theta\mathbf{e}_{y}+\Omega\sin\theta\mathbf{e}_{z}$

with respect to the rotating reference frame. Moreover, given any rotating
orthonormal basis $\\{\mathbf{e}_{x},\mathbf{e}_{y},\mathbf{e}_{z}\\}$ with
rotation vector $\mathbf{\Omega}$ it is always the case that the basis vectors
satisfy

(3)
$\displaystyle\frac{d}{dt}{\mathbf{e}_{k}}=\mathbf{\Omega}\times\mathbf{e}_{k},\quad\text{for
}k\in\\{x,y,z\\}.$

and the reader is referred to [20] for a comprehensive derivation of this
formula. Thus a fluid particle following some trajectory

${\bf r}(t)=X(t)\mathbf{e}_{x}+Y(t)\mathbf{e}_{y}+Z(t)\mathbf{e}_{z},$

will have a corresponding velocity and acceleration given by

(4) $\mathbf{u}=\frac{\partial\mathbf{r}}{\partial
t}+\mathbf{\Omega}\times\mathbf{r}\qquad\mathbf{a}=\frac{\partial\mathbf{u}}{\partial
t}+2\mathbf{\Omega}\times\mathbf{u}+\mathbf{\Omega}\times\mathbf{\Omega}\times\mathbf{r}.$

The acceleration of a fluid element, $\mathbf{a}$, has two contributions due
to fictitious forces arising from the rotation of the Earth. The terms
proportional to $\Omega^{2}$ correspond to the centripetal acceleration, which
is the acceleration of a body away from the surface of the Earth due to the
rotation of the surface. In practice, this term is negligible since it is
effectively cancelled by the gravitational acceleration of the object,
although, the effects of this force are observable in the deformation of the
Earth away from a perfect sphere (see [16] for further discussion). The term
proportional to $\Omega$ is the result of the Coriolis force, and in the
following we shall only consider the effect of this fictitious force.

####  2.1.1. The Euler equation in rotating coordinates

Given equation (4) and the effective cancellation of centripetal and
gravitational forces on the surface of the Earth, it is found that the
familiar material derivative is modified according to

$\frac{\partial}{\partial t}+\mathbf{u}\cdot\nabla\to\frac{\partial}{\partial
t}+\mathbf{u}\cdot\nabla+2\mathbf{\Omega}\times,$

thus modifying the Euler equation according to

(5) $\begin{cases}u_{t}+uu_{x}+vu_{y}+wu_{z}+2\Omega\cos\theta
w-2\Omega\sin\theta v=-\frac{1}{\rho}P_{x}\\\
v_{t}+uv_{x}+vv_{y}+wv_{z}+2\Omega\sin\theta u=-\frac{1}{\rho}P_{y}\\\
w_{t}+uw_{x}+vw_{y}+ww_{z}-2\Omega\cos\theta
u=-\frac{1}{\rho}P_{z}-g.\end{cases}$

The hydrodynamic pressure in the fluid is represented by $P$ with the fluid
density denoted by $\rho$, while $u$, $v$ and $w$ denote the velocity field
components.

In the following, the meridional coordinate is restricted to the interval
$-5^{\circ}<\theta<5^{\circ}$, and in this equatorial region the approximation

$\cos\theta\simeq 1\qquad\sin\theta\simeq\frac{y}{R}$

may be used in equation (5), cf. [7, 10, 11, 13]. Introducing the parameter

$\beta=\frac{2\Omega}{R}$

this approximation yields the so-called $\beta$-plane formulation, given by

(6) $\begin{cases}u_{t}+uu_{x}+wu_{z}+2\Omega w=-\frac{1}{\rho}P_{x}\\\ \beta
yu=-\frac{1}{\rho}P_{y}\\\ w_{t}+uw_{x}+ww_{z}-2\Omega
u=-\frac{1}{\rho}P_{z}-g,\end{cases}$

where the velocity field is restricted by $v=0$ and $u_{y}=w_{y}=0$, meaning
the equatorial flows under consideration are effectively restricted to zonal
motion. While there is no motion in the meridional direction, this does not
exclude variation in that direction, and in particular the depth of each
interface between the layers of the flow will depend on latitude, thus making
the flow three-dimensional. Incompressibility of the flow is ensured when

(7) $u_{x}+v_{y}+w_{z}=0,$

while mass-conservation along the flow follows from

(8) $\rho_{t}+\nabla\cdot(\mathbf{u}\rho)=0,$

and the system (6)–(8) constitute the governing equations of the model. The
boundary conditions coupled to these governing equations are of the form

(9) $\displaystyle\begin{rcases}P=P_{0}\\\
w=\eta_{t}+u\eta_{x}\end{rcases}\quad\text{on }z=\eta(x,t)$
$\displaystyle(u,w)=(0,0)\quad\text{for
}z<-D+\frac{\kappa\beta}{4\Omega}y^{2},$

where $\kappa>0$ is constant, and whose relevance will be clarified in section
4. The first of these boundary conditions ensures the pressure on the free
surface $z=\eta(x,t)$ is given by the constant atmospheric pressure $P_{0}$.
The second boundary condition on this free surface implies that a fluid
particle always remains on this surface. The final kinematic boundary
condition ensures the fluid motion ceases below a certain depth, which
demarcates the motionless layer.

##  3\. The flow about the equator

The model resembles that presented in [11], in that the waves occur in a fluid
with two regions of different densities. The interface where the density
changes is the thermocline, and it is this change in density which allows
internal waves to arise. These internal waves propagate towards the east and
move with constant wave-speed $c$ which will be prescribed by the dispersion
relation in equation (30). The waves, as already pointed out, have zero
meridional velocity, and occur beneath a surface layer $\mathcal{L}(t)$ whose
behaviour is heavily influenced by wind effects. The behaviour of the fluid in
this layer is not under consideration in this work. The lower boundary of the
layer is described by the surface $z=\eta_{+}(x-ct,y)$, and beneath this
boundary there are four regions wherein the fluid motion exhibits distinctive
behaviour.

The upper most layer we investigate in this work is given by

(10)
$\mathcal{M}(t)=\left\\{(x,y,z):x\in\mathbb{R},y\in\mathbb{R},z\in\left(\eta_{0}(x-ct,y),\eta_{+}(x-ct,y)\right)\right\\}.$

Wind effects on the internal flow in this work will be neglected for
simplicity, however, it is known that atmospheric motion can play a role in
the dynamics of internal flows in stratified layers, cf. [13, 32]. The fluid
density within this region is of constant value $\rho_{0}$, while the lower
boundary of this region corresponds to the thermocline. Beneath the
thermocline, the uniform flow layer is bounded from below by the surface

$z=-d+\frac{\kappa\beta}{4\Omega}y^{2}$

where $d$ is constant. In this region the fluid density increases slightly to
$\rho_{+}>\rho_{0}$. The typical relative change in fluid density observed in
the equatorial Pacific is of the order

(11) $\frac{\rho_{+}-\rho_{0}}{\rho_{0}}\approx 4\times 10^{-3},$

(see [37]). Beneath this uniform flow layer, the transitional layer is bounded
below by the surface

$z=-D+\frac{\kappa\beta}{4\Omega}y^{2},$

wherein the fluid is again of constant uniform density $\rho_{+}$. The final
layer, is unbounded below and corresponds to the motionless layer, wherein the
fluid motion vanishes completely, and the fluid density is again $\rho_{+}$
throughout.

Figure 2. A cross section of the flow in the equatorial plane $y=0$ with a
flow wavelength of $L=200\,\mathrm{m}$. The average depth of the near surface
layer is $60\,\mathrm{m}$, while the thermocline lies at an average depth of
$120\,\mathrm{m}$. The transitional layer begins at $160\,\mathrm{m}$
approximately, while the motionless layer begins at $200\,\mathrm{m}$ beneath
the free surface.

The layers of the fluid domain and the boundaries between those layers are
shown in Figure 2. In particular, the thermocline exhibits a cusp at the wave
trough, which is one of the distinctive features of extreme waves. In the case
of extreme waves, the curve described by the thermocline is a cycloid, which
is the path traced out by a particle at a fixed point on the circumference of
a circle which rolls without slipping. In the case of smooth internal waves
investigated in [11], the thermocline wave profile is a trochoid, which is the
path traced out by some point along the radius of a circle which rolls without
slipping, cf. [6].

###  3.1. The layers beneath the thermocline

The behaviour of the flow in the layers beneath the thermocline is identical
to that in [11], and we outline here the relevant features of these layers for
the sake of completeness. The significant difference between the wave of
extreme form and the internal waves previously investigated in that work occur
in the layer $\mathcal{M}(t)$, in particular along the thermocline, due to the
presence of cusps at the wave-trough as illustrated in Figure 2.

####  3.1.1. The motionless layer

At any meridional location $y$ the upper boundary of the motionless layer is
given by

$z=-D+\frac{\kappa\beta}{4\Omega}y^{2},$

and beneath this boundary the fluid is still with $u=v=w=0$. In this case, the
hydro-static pressure is given by

(12) $P=P_{0}-\rho_{+}gz\quad\text{when
}z\leq-D+\frac{\kappa\beta}{4\Omega}y^{2}$

The constant parameter $D$ is the depth of the upper boundary of this layer,
as measured along the equator $y=0$.

####  3.1.2. The transitional layer

The upper boundary of this region is given by the surface

$z=-d+\frac{\kappa\beta}{4\Omega}y^{2},$

where the constant parameter $d$ is the depth of the boundary as measured
along the equator. Within this region the fluid is motionless in the vertical
and meridional directions, however, the fluid does move in the zonal
direction, with a velocity which depends on the depth and the meridional
location. Specifically, the horizontal velocity in this layer is given by

(13) $\displaystyle
u(x,y,z)=\frac{c}{D-d}\left(z-\frac{\kappa\beta}{4\Omega}y^{2}+D\right),\quad
v=0,\quad w=0.$

The vertical component of the Euler equation within this transitional layer
now becomes

(14) $\frac{2\Omega
c\rho_{+}}{D-d}\left(z-\frac{\kappa\beta}{4\Omega}y^{2}+D\right)=P_{z}+\rho_{+}g,$

and upon integrating we have

(15) $\displaystyle P$
$\displaystyle=P|_{z=-D+\frac{\kappa\beta}{4\Omega}y^{2}}+\left[-\rho_{+}gz^{\prime}+\frac{2\Omega
c\rho_{+}}{D-d}\left(\frac{z^{\prime
2}}{2}-\frac{\kappa\beta}{4\Omega}z^{\prime}y^{2}+Dz^{\prime}\right)\right]^{z}_{-D+\frac{\kappa\beta}{4\Omega}y^{2}}$
$\displaystyle=P_{0}-\rho_{+}gz+\frac{\Omega
c\rho_{+}}{D-d}\left(z+D-\frac{\kappa\beta}{4\Omega}y^{2}\right)^{2}.$

The pressure on the lower boundary of this region,
$P|_{z=-D+\frac{\kappa\beta}{4\Omega}y^{2}}$, is deduced from the continuity
of the pressure across the interface of two layers in the fluid domain, along
with the expression (12) for the pressure in the motionless layer.

####  3.1.3. The uniform flow layer

The upper boundary of the uniform flow layer corresponds to the thermocline

$z=\eta_{0}(x-ct,y),$

where $\eta_{0}$ is an eastward moving wave-form, moving with constant speed
$c$ along any fixed meridional line $y$. Within the layer

$-d+\frac{\kappa\beta}{4\Omega}y^{2}\leq z<\eta_{0}(x-ct,y)$

the flow is uniform, given by

(16) $u=c,\quad v=0,\quad w=0.$

Replacing this velocity field in the Euler equation (6), we deduce that

(17) $P(x,y,z)=P_{0}-\rho_{+}gz+2\rho_{+}\Omega
c\left(z+\frac{D+d}{2}-\frac{\kappa\beta}{4\Omega}y^{2}\right),$

having integrated each component and imposed equation (15) on the lower
boundary $z=-d+\frac{\kappa\beta}{4\Omega}y^{2}$.

###  3.2. The layer above the thermocline

The layer immediately above the thermocline, given by

$\eta_{0}(x-ct,y)\leq z\leq\eta_{+}(x-ct,y)$

is fully dynamic, in that the velocity field does not reduce to a simplified
form as in all the layers beneath the thermocline. In this layer, it is
convenient to analyse the flow in terms of Lagrangian coordinates, $(q,s,r)$,
which serve to label the individual fluid particles within this layer. The
flow in this layer is described in parametric form by

(18)
$\begin{cases}X(t;q,s,r)=q-\frac{1}{k}e^{-k\left(r+f(s)\right)}\sin(k(q-ct))\\\
Y(t;q,s,r)=s\\\
Z(t;q,s,r)=-d_{0}+r-\frac{1}{k}e^{-k(r+f(s))}\cos(k(q-ct)),\end{cases}$

where $L=\frac{2\pi}{k}$ is the wavelength of the flow, while the wave speed
$c$ is deduced from the continuity of the pressure across the thermocline.

Figure 3. The circular paths in the equatorial plane $y=0$ traced by particles
in the layer $\mathcal{M}(t)$ are circles with centre $(q,r-d_{0})$ whose
radius increases with depth. The wavelength of the solution shown is $L=200\
\mathrm{m}$. The particle trajectories are counterclockwise and the particle
positions shown in the figure are at time $t=0.$

The function $f:=f(s)$ controls the decay of the fluid disturbance away from
the equator which corresponds to the line $s=0$. The parameter space variables

$(q,s,r)\in\mathbb{R}\times(-s_{0},s_{0})\times(r_{0}(s),r_{+}(s))$

may be interpreted as the initial locations of a motionless fluid body prior
to being disturbed by a remotely generated wave propagating in the form (18).
The functions $r_{0}(s)$ and $r_{+}(s)$ characterise the surfaces
$z=\eta_{0}(x,y,t)$ and $z=\eta_{+}(x,y,t)$ respectively, while $s_{0}\leq
250\,\mathrm{km}$ is the typical meridional extent of such disturbances (cf.
[11]). Under the action of this flow, the particle labelled $(q,s,r)$
describes a circular path in the plane of fixed $s$, with centre $(q,r-d_{0})$
and radius $\frac{1}{k}e^{-k(r+f(s))}$ which decreases with depth, as
illustrated in Figure 3.

The dynamic possibility of the flow (18) has been confirmed in the context of
surface waves in the work [38], and in the context of internal geophysical
waves in [35]. In either regime, as time advances, the mapping (18) remains a
diffeomorphism from the parameter space $(q,s,r)$ to the fluid domain when
$r\geq 0$. Meanwhile, the wave profile described by (18), at any fixed
$(s,r,t)$ is characterised to a large extent by the radius $e^{-k(r+f(s))}$.
Crucially, for this profile to correspond to the graph of a function it is
necessary that $r+f(s)\geq 0$, with the wave of extreme form emerging when
$r+f(s)=0.$ When $r+f(s)>0$ the wave-form is trochoidal, while the case
$r+f(s)<0$ yields a curve which intersects with itself and so does not
correspond to the graph of a function, cf. [6, 11].

Introducing the parameters

$\xi=k(r+f(s))\quad\theta=k(q-ct),$

the Jacobian matrix of the coordinate transform (18) is of the form

(19)
$\frac{\partial(X,Y,Z)}{\partial(q,s,r)}=\begin{bmatrix}1-e^{-\xi}\cos\theta&0&e^{-\xi}\sin\theta\\\
f_{s}e^{-\xi}\sin\theta&1&f_{s}e^{-\xi}\cos\theta\\\
e^{-\xi}\sin\theta&0&1+e^{-\xi}\cos\theta\end{bmatrix}$

while corresponding Jacobian is given by

(20) $J=1-e^{-2\xi}.$

We note that this Jacobian is time-independent indicating that the flow is
volume preserving. In the case of smooth waves studied in [11], the Jacobian
satisfies $J>0$ at all points in the closure of the layer $\mathcal{M}(t)$,
and this allows the use of the implicit function theorem to infer the
existence of an associated profile $r_{0}(s)$. In contrast, for extreme waves
this strict inequality now becomes $J\geq 0$, with the Jacobian vanishing on
the thermocline and specifically along the equator, which gives rise to a
number of mathematical complications. In particular it precludes the use of
the implicit function theorem to confirm the existence of waves of extreme
form propagating in the thermocline, an issue which is addressed in the
Section 3.3.

A particularly appealing feature of the Lagrangian formalism is that the fluid
velocity at any point along the flow given in equation (18) follows directly
from

(21) $\begin{cases}u=\frac{DX}{Dt}=ce^{-\xi}\cos\theta\\\ v=\frac{DY}{Dt}=0\\\
w=\frac{DZ}{Dt}=-ce^{-\xi}\sin\theta.\end{cases}$

The corresponding acceleration along the flow is similarly obtained, yielding

(22) $\begin{cases}\frac{Du}{Dt}=kc^{2}e^{-\xi}\sin\theta\\\
\frac{Dv}{Dt}=0\\\ \frac{Dw}{Dt}=kc^{2}e^{-\xi}\cos\theta.\end{cases}$

In terms of the labelling variables $(q,s,r)$, the gradient of the
hydrodynamic pressure $P$ is obtained via

(23)
$\nabla_{(q,s,r)}P=\frac{\partial(X,Y,Z)}{\partial(q,s,r)}\nabla_{(x,y,z)}P,$

and so applying equations (19)–(23) to the $\beta$-plane approximation (6), we
find

(24) $\begin{cases}P_{q}=-\rho_{0}(kc^{2}-2\Omega c+g)e^{-\xi}\sin\theta\\\
P_{s}=-\rho_{0}(kc^{2}-2\Omega c)f_{s}e^{-2\xi}-\rho_{0}(\beta
cs+f_{s}g)e^{-\xi}\cos\theta\\\ P_{r}=-\rho_{0}(kc^{2}-2\Omega
c)e^{-2\xi}-\rho_{0}(kc^{2}-2\Omega
c+g)e^{-\xi}\cos\theta-\rho_{0}g.\end{cases}$

Integrating the expression for $P_{q}$ and comparing with the expression for
$P_{s}$, we deduce that

(25) $f(s)=\frac{\beta s^{2}}{2(kc-2\Omega)}.$

Moreover the pressure distribution, written in terms of the labelling
variables, is given by

(26) $P(q-ct,s,r)=\rho_{0}\frac{kc^{2}-2\Omega
c}{2k}e^{-2\xi}-\rho_{0}gr+\rho_{0}\frac{kc^{2}-2\Omega
c+g}{k}e^{-\xi}\cos\theta+\tilde{P}_{0},$

where $\tilde{P}_{0}$ is a constant.

####  3.2.1. The wave phase speed

Since the thermocline $z=\eta_{0}(x,y,t)$ is identified with the parameterised
surface $r=r_{0}(s)$, the pressure in the uniform current layer given by
equation (17) combined with the third component of equation (18) evaluated on
$r_{0}(s)$, give the following expression for the pressure on the thermocline

(27) $P(q-ct,s,\eta_{0})=P_{0}-\rho_{+}\frac{\kappa\beta
c}{2}s^{2}+\rho_{+}(2\Omega
c-g)\left(-d_{0}+r_{0}(s)-\frac{1}{k}e^{-\xi_{0}}\cos(\theta)\right)$

where $\xi_{0}=k(r_{0}(s)+f(s))$. Alternatively, the pressure in the layer
$\mathcal{M}(t)$ given by (26) and restricted to the thermocline $r=r_{0}(s)$
also yields

(28) $\displaystyle P(q-ct,s,r_{0}(s))$
$\displaystyle=\tilde{P}_{0}+\rho_{0}\frac{kc^{2}-2\Omega
c+g}{k}e^{-\xi_{0}}\cos\theta$
$\displaystyle\qquad+\rho_{0}\frac{kc^{2}-2\Omega
c}{2k_{0}}e^{-2\xi}-\rho_{0}gr_{0}(s).$

Since the pressure is continuous throughout the flow at all times, it follows
from expressions (27)–(28) that

(29) $\rho_{+}(g-2\Omega c)=\rho_{0}(kc^{2}-2\Omega c+g),$

from which we conclude

(30)
$c=\frac{\rho_{+}-\rho_{0}}{k\rho_{0}}\left(\sqrt{\Omega^{2}+kg\frac{\rho_{0}}{\rho_{+}-\rho_{0}}}-\Omega\right)>0,$

having retained only the positive root since the waves travel eastward.
Crucially this dispersion relation highlights the necessity of a heterogeneous
fluid body for the existence of these internal waves, since it is immediately
clear that if $\rho_{+}=\rho_{0}$, then $c=0$ and the internal waves no longer
exist. Decay of the internal wave with meridional distance from the equator is
assured when $kc-2\Omega>0$, and so it follows from (30) that

(31)
$k>\frac{4\Omega^{2}}{g}\left(1+\frac{\rho_{0}}{\rho_{+}-\rho_{0}}\right),$

which combined with inequality (11), ensures that a maximal wavelength for
internal waves is of order $L=\frac{2\pi}{k}\approx 1.2\times
10^{7}\mathrm{m}.$

###  3.3. Existence of extreme waves in the thermocline

Equating the expressions given in equations (27)–(28) ensures the continuity
of the pressure across the interface of the the layer $\mathcal{M}(t)$ and the
uniform current layer. The equivalence of these two expression allows us to
deduce an implicit relation for the function $r_{0}(s)$ of the form

(32) $\rho_{+}\frac{\kappa\beta c}{2}s^{2}+\rho_{0}(kc^{2}-2\Omega
c)\left(\frac{1}{2k}e^{-2k(r_{0}(s)+f(s))}+r_{0}(s)\right)\\\
=P_{0}-\tilde{P}_{0}+\rho_{+}gd_{0}-\rho_{+}\Omega c(2d_{0}-D-d).$

Let us denote

$G(s.r)=\rho_{+}\frac{\kappa\beta c}{2}s^{2}+\rho_{0}(kc^{2}-2\Omega
c)\left(\frac{1}{2k}e^{-2k(r+f(s))}+r\right),$

and differentiating with respect to $r$ we find

(33) $G_{r}(s,r)=\rho_{0}(kc^{2}-2\Omega c)(1-e^{-2k(r+f(s))})>0,\text{ when
}r>0.$

In particular, smooth nonlinear flows always satisfy $r_{0}(s)\geq
r_{0}^{*}>0$ for all $s\in(-s_{0},s_{0})$, and so under these circumstances it
is clear that $r\mapsto G(s,r)$ is a strictly increasing diffeomorphism from
the open interval $(0,\infty)$ onto the open interval
$\left(\rho_{+}\frac{\kappa\beta c}{2}s^{2}+\rho_{0}\frac{(kc^{2}-2\Omega
c)}{2k}e^{-2kf(s)},\infty\right)$. This in turn ensures the implicit function
theorem may be used to infer the existence of a solution $r=r_{0}(s)$ of
equation (32). In addition we observer that $G(s,r)$ is effectively a function
of $s^{2}$ only, and so the implicit function theorem ensures $r_{0}(s)$ is
effectively a function of $s^{2}$ also, in which case it is clear that
$r_{0}(s)$ is an even function of $s$, see [11] for further discussion.

Extreme waves emerge from (18) only when $r_{0}(s)\geq r_{0}^{*}=0$ with
$r_{0}(s)=0$ along the equator $s=0$ (in contrast $r_{0}^{*}>0$ for smooth
waves). In this scenario equation (33) now becomes $G_{r}(s,r)\geq 0$ for
$r\geq r_{0}$ and so significantly $G(s,r)$ is no longer a strictly increasing
function of $r$ for all $s\in(-s_{0},s_{0})$ (see also the discussion after
equation (20)). Equivalently, the fact that $G_{r}(s,r)$ is not strictly
positive everywhere essentially means $G(s,r)$ is not automatically an
injective function of $r$ for each fixed $s\in(-s_{0},s_{0})$, and so the
implicit function theorem cannot be automatically applied to infer a solution
of the relation (32). Nevertheless, even in the case of extreme waves it can
still be shown for any fixed $s\in(-s_{0},s_{0})$ that the function $r\mapsto
G(s,r)$ is an injective map as follows:

######  Lemma 3.1.

Given any fixed $s$ the function

$\displaystyle G(s,\cdot):[0,\infty)\to\left[\rho_{+}\frac{\kappa\beta
c}{2}s^{2}+\rho_{0}\frac{(kc^{2}-2\Omega c)}{2k}e^{-2kf(s)},\infty\right)$
$\displaystyle G(s,\cdot):r\mapsto\rho_{+}\frac{\kappa\beta
c}{2}s^{2}+\rho_{0}(kc^{2}-2\Omega
c)\left(\frac{1}{2k}e^{-2k(r+f(s))}+r\right)$

is injective.

###### Proof.

Away from the equator we see that

$G_{r}(s,r)=\rho_{0}(kc^{2}-2\Omega c)(1-e^{-2k(r+f)})>0,$

since $\frac{\beta s^{2}}{2(kc-2\Omega)}>0$ when $s\neq 0$ and hence $r\mapsto
G(s,r)$ is a strictly increasing function when $r\in[0,\infty)$.

When $s=0$ we note that $G(0,r)$ is continuous for any $r\in[0,\infty)$ and
differentiable for any $r\in(0,\infty)$. Hence, the mean value theorem ensures
that for any pair of values $0\leq r_{1}<r_{2}<\infty$, there exists some
$\gamma\in(r_{1},r_{2})$ such that

(34) $\displaystyle\frac{G(0,r_{2})-G(0,r_{1})}{r_{2}-r_{1}}$
$\displaystyle=\rho_{0}(kc^{2}-2\Omega c)\left(1-e^{-2k(\gamma)}\right)>0$

since $0\leq r_{1}<\gamma$ and $kc^{2}-2\Omega c>0$. Thus $G(0,r)$ is
injective for $r\in[0,\infty).$ It follows that the function $G(s,r)$ is an
injective function of $r$ for every fixed $s$. ∎

Using the fact that $G(s,r)$ is an injective function on the interval
$r\in[0,\infty)$ we may now show that there exists a solution of the implicit
relation (32).

######  Theorem 3.2.

There exists an $s_{0}>0$ such that for any $s\in(-s_{0},s_{0})$ the implicit
relation

$G(s,r)=P_{0}-\tilde{P}_{0}+\rho_{+}gd_{0}-\rho_{+}\Omega c(2d_{0}-D-d)$

has a unique solution $r=r_{0}(s)\geq 0$ such that $r(0)=0$. Moreover,
$r_{0}(s)$ is a continuous function for $s\in(-s_{0},s_{0})$.

###### Proof.

Firstly, we note since the right-hand side of equation (32) is constant and
that $r_{0}(0)=0$, then equivalently we want to show there exists a continuous
solution $r_{0}(s)$ to the implicit relation $G(s,r)=G(0,0)$. Let us consider
the behaviour of $G$ along the line segment

$\left\\{(0,r):0\leq r<\infty\right\\}.$

Figure 4. A pictorial outline for the proof of Theorem 3.2, confirming the
existence of a solution of the implicit relation $G(s,r)=G(0,0)$ of extreme
form.

It follows from Lemma 3.1 that for any $r_{1}>0$ then $G(0,r_{1})>G(0,0)$.
Noting that $G(s,r)$ is continuous everywhere, it follows that an interval
$(-s_{0},s_{0})$ may always be chosen so that

(35) $G(s,r_{1})>G(0,0),\text{ for }s\in(-s_{0},s_{0}),$

see Figure 4. Now fixing some $s_{1}\in(-s_{0},s_{0})$ with $s_{1}\neq 0$, we
have

(36) $G(s_{1},r_{1})>G(0,0)$

again. Moreover, since $G(s_{1},r)$ is a continuous and injective function of
$r$ (cf. Lemma 3.1), it follows that there is a unique
$r_{0}(s_{1})\in(0,r_{1})$ such that

(37) $G(s_{1},r_{0}(s_{1}))=G(0,0),$

as illustrated in Figure 4. Hence, for all $s\in(-s_{0},s_{0})$, there exists
a function $r_{0}(s)$ such that $G(s,r_{0}(s))=G(0,0)$.

The continuity of $r_{0}(s)$ may be seen to follow in a similar way. Fixing
some $s_{1}\in(-s_{0},s_{0})$, to confirm the continuity we must show that for
an arbitrarily small $\epsilon$ there is a a corresponding $\delta>0$ such
that

$\left|{r_{0}(s)-r_{0}(s_{1})}\right|<\epsilon,\text{ for all
}s\in(s_{1}-\delta,s_{1}+\delta).$

However, we have already shown that an interval $(-s_{0},s_{0})$ may always be
chosen which ensures the existence of some function $r(s)$ for each
$s\in(-s_{0},s_{0})$. Consequently we may always choose an appropriate sub-
interval $(s_{1}-\delta,s_{1}+\delta)\subset(-s_{0},s_{0})$ for each
$s_{1}\in(-s_{0},s_{0})$, so that the continuity of $r(s)$ is ensured (see
[15] for further discussion). ∎

##  4\. Physical features of the flow

###  4.1. The profile near the equator

Given any point $(s,r_{0}(s))$ where $G_{r}(s,r_{0}(s))\break>0$, the
derivative of $r_{0}(s)$ with respect to $s$ is given by

(38)
$r^{\prime}(s)=-\left.\frac{G_{s}(s,r)}{G_{r}(s,r)}\right|_{r=r_{0}(s)}=-\frac{\beta
cs(\rho_{+}\kappa-\rho_{0}e^{-2\xi_{0}})}{\rho_{0}(kc^{2}-2\Omega
c)(1-e^{-2\xi_{0}})}.$

We recall that Lemma 3.1 and Theorem 3.2 require $r\in[0,\infty)$, and so it
follows that $r_{0}(s)$ must not decrease, at least initially, as we move away
from the equator. Hence, it follows from equation (38) that
$r_{0}^{\prime}(s)\leq 0$ for $s<0$ and $r_{0}^{\prime}(s)\geq 0$ for $s>0$,
at least about the equator. This in turn requires

$\rho_{+}\kappa-\rho_{0}e^{-2\xi_{0}}\leq 0,\text{ for all }s\neq 0$

from which it follows that

(39) $\kappa<\frac{\rho_{0}}{\rho_{+}}<1.$

We note from equation (38) if $\kappa=\frac{\rho_{0}}{\rho_{+}}$ then
$r_{0}^{\prime}(s)\leq 0$ when $s\geq 0$ and vice-versa, meaning $r_{0}(s)$
decreases in moving away from the equator, hence the necessity for a strict
inequality in (39).

When $\kappa<\frac{\rho_{0}}{\rho_{+}}$, equation (32) has an exact solution
of the form

(40)
$r(s)=r_{0}^{*}+\frac{1}{2k}-\kappa\frac{\rho_{+}}{\rho_{0}}\frac{\beta}{(kc-2\Omega)}s^{2}+W\left(-e^{-1-2k(1-\kappa\frac{\rho_{+}}{\rho_{0}})\frac{\beta}{kc-2\Omega}s^{2}}\right),$

where $W(\cdot)$ is the Lambert $W$ function, defined by

$z=W(z)e^{W(z)}\quad\text{ for all }z\in\mathbb{C},$

The function $W(x)$ is real valued for any $x\in[-\frac{1}{e},\infty)$. The
function is double valued on the interval $[-\frac{1}{e},0)$ with an
increasing and decreasing branch in this interval, while the function is
single valued on the interval $[0,\infty).$ An interesting application of the
Lambert $W$ function for a convective Rayleigh-Benard model may be found in
[41]. In Figure 5 in the left panel, the graph of $(s,r_{0}(s))$ is shown for
a wave length $L=200\,$m and as expected the function is initially increasing
when moving away from the equator. In the right panel, the cross section
$(Y(q,s,r_{(}s),t),Z(q,s,r_{0}(s),t))$ is illustrated and again we see the
shape of $r_{0}(s)$ reflected in this graph, except at the wave crest where
the profile becomes smooth. In Figure 6 the thermocline is illustrated from
beneath to emphasise the cusp at the wave trough, one of the hallmarks of
extreme waves. We also note from the figures that the wave profile does not
appear to be differentiable at the wave trough, again a common feature of
extreme waves.

Figure 5. The left-hand panel shows the graph $(s,r_{0}(s))$ when $\kappa=0.5$
and $L=200$ meters. The parameter $d_{0}=120$ ensures a mean equatorial depth
of $120\,$m. In the right-hand panel the corresponding thermocline profile in
the fluid domain, at the equatorial locations indicated in the figure. Figure
6. The thermocline surface when $L=200\,$m, $r_{0}^{*}=0\,$m and
$d_{0}=120\,$m, with $\kappa=0.5$. The surface is viewed from beneath and the
value of $\kappa$ has been chosen to emphasise the extreme behaviour of the
wave along the equator at the wave-troughs.

When the wavelength $L$ is sufficiently small, there are values
$s\in(-s_{0},s_{0})$ such that the solution $r_{0}(s)<0$, for instance when
$L=100\,$m we find that $r_{0}(s)<0$ when $|s|\gtrsim 234\,$km. However, we
also note that for $|s|\gtrsim 234\,$km then $f(s)+r(s)>0$, meaning the
details of Lemma 3.1 and Theorem 3.2 remain valid, and crucially the graph of
$(X(q,r_{0}(s),s,t),Z(q,r_{0}(s),s,t))$ does not intersect itself. Thus, it is
only near the equator that we must insist $r_{0}(s)>0$, in line with our
considerations above.

###  4.2. Current flow reversal

In [11] it was shown that the mean Eulerian velocity of the solution (18)
induces a westward flowing current in the near surface region $M(t)$, and
owing to the fact that current in the layer beneath the thermocline is
eastward, it was found that such behaviour reflects the features of the
equatorial undercurrent (with relevant data in [42]). We shall now demonstrate
that this property of the dynamic solution is also a feature of the flow in
the extreme case, thus maintaining one of the most attractive features of the
flow from a geophysical perspective.

At any given latitude $s\in(-s_{0},s_{0})$, any fixed depth $z=z_{0}$ within
the layer $M(t)$ may be written in terms of Lagrangian variables via the
solution (18), in the form

(41) $z_{0}=-d_{0}+r+\frac{1}{k}e^{-k(r+f(s))}\cos(k(q-ct)),$

whose solution is defined implicitly via

(42) $r=\zeta(q,s,z_{0},t).$

Since the thermocline $r(s)=r_{0}^{*}$ is the lower boundary of this layer, it
follows that $\zeta(q,s,z_{0},t)\geq r_{0}^{*}$ for all $q,s,t$ and
$z\geq\eta_{0}(x,y,t)$. Replacing (42) in (41) and differentiating with
respect to $q$ we deduce

(43) $\zeta_{q}=-\frac{e^{-\xi}\sin\theta}{1+e^{-\xi}\cos\theta}.$

The mean horizontal flow at this latitude $s$ and depth $z_{0}$ is then given
by

(44)
$c+\bar{u}(s,z_{0})=\frac{1}{T}\int_{0}^{T}\left[c+u(x-ct,s,z_{0})\right]dt.$

Using

$t=\frac{L}{c}\qquad dx=cdt$

along with the change of variables (having fixed $s$ and $z$)

$dx=\frac{\partial x}{\partial q}dq,$

it follows from (18), (21) and (43) that the mean horizontal current is given
by the integral

(45)
$\bar{u}(s,z_{0})=-\frac{c}{L}\int_{0}^{L}e^{-2\xi}dq\Rightarrow-c<\bar{u}(s,z_{0})<0,$

meaning there is a net westward flow in the layer $M(t)$. Thus there is a flow
reversal in the transition across the thermocline when the profile is smooth
or of extreme form.

## Acknowledgments

The author is grateful to the referees for several helpful comments. The
author is grateful to the organisers of the programme “Mathematical Aspects of
Physical Oceanography (January 2019 – March 2019)” at the Erwin Schrödinger
Institute, Vienna, Austria.

## References

  * [1] C. J. Amick, L. E. Fraenkel and J. F. Toland,  On the Stokes conjecture for the wave of extreme form,  _Acta Math._ , 148 (1982), 193–214.
  * [2] A. Bennett,  _Lagrangian Fluid Dynamics_ ,  Cambridge Monographs on Mechanics, Cambridge University Press, 2006.
  * [3] B. Buffoni and J. Toland,  _Analytic Theory of Global Bifurcation_ ,  Princeton University Press, 2003.
  * [4] A. Compelli and R. Ivanov,  On the dynamics of internal waves interacting with the equatorial undercurrent,  _J. Nonlin. Math. Phys._ , 22 (2015), 531–539.
  * [5] A. Constantin,  On the deep water wave motion,  _J. Phys. A_ , 34 (2001), 1405–1417.
  * [6] A. Constantin,  _Nonlinear Water Waves with Applications to Wave-Current Interactions and Tsunamis_ , vol. 81 of CBMS-NSF Regional Conference Series in Applied Mathematics,  SIAM, Philadelphia, 2011.
  * [7] A. Constantin,  An exact solution for equatorially trapped waves,  _J. Geophys. Res. Oceans_ , 117 (2012).
  * [8] A. Constantin,  On the modelling of equatorial waves,  _Geophys. Res. Lett._ , 39 (2012).
  * [9] A. Constantin,  Particle trajectories in extreme Stokes waves,  _IMA J. Appl. Math._ , 77 (2012), 293–307.
  * [10] A. Constantin,  Some three-dimensional nonlinear equatorial flows,  _J. Phys. Oceanog._ , 43 (2013), 165–175.
  * [11] A. Constantin,  Some nonlinear, equatorially trapped, nonhydrostatic internal geophysical waves,  _J. Phys. Oceanogr._ , 44 (2014), 781–789.
  * [12] A. Constantin and P. Germain,  Instability of some equatorially trapped waves,  _J. Geophys. Res. Oceans_ , 118 (2013), 2802–2810.
  * [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 S. G. Monismith,  Gerstner waves in the presence of mean currents and rotation,  _J. Fluid Mech._ , 820 (2017), 511–528.
  * [15] R. Courant,  _Differential and Integral Calculus_ ,  John Wiley & Sons, Inc., New York, 1988.
  * [16] B. Cushman-Roisin and J. M. Beckers,  _Introduction to Geophysical Fluid Dynamics: Physical and Numerical Aspects_ , vol. 101,  Academic Press, 2011.
  * [17] M. L. Dubreil-Jacotin,  Sur les ondes de type permanent dans les liquides,  _Atti. Accad. Naz. Lincei, Mem. Cl. Sci. Fis. Mat. Nat._ , 15 (1932), 814–819.
  * [18] F. Genoud and D. Henry,  Instability of equatorial water waves with an underlying current,  _J. Math. Fluid Mech._ , 16 (2014), 661–667.
  * [19] F. Gerstner,  Theorie der wellen samt einer daraus abgeleiteten theorie der deichprofile,  _Ann. Phys._ , 2 (1809), 412–445.
  * [20] H. Goldstein, C. P. Poole and J. L. Safko,  _Classical Mechanics_ ,  Pearson International Edition, Addison Wesley, 2002.
  * [21] D. Henry,  On Gerstner’s water wave,  _J. Nonlin. Math. Phys._ , 15 (2008), 87–95.
  * [22] D. Henry,  An exact solution for equatorial geophysical water waves with an underlying current,  _Eur. J. Mech. B Fluids_ , 38 (2013), 18–21.
  * [23] D. Henry,  On three-dimensional Gerstner-like equatorial water waves,  _Phil. Trans. R. Soc. A_ , 376 (2018), 20170088, 16pp.
  * [24] D. Henry and H. C. Hsu,  Instability of internal equatorial water waves,  _J. Diff. Eqn._ , 258 (2015), 1015–1024.
  * [25] D. Henry and S. Sastre-Gomez,  Mean flow velocities and mass transport for equatorially-trapped water waves with an underlying current,  _J. Math. Fluid Mech._ , 18 (2016), 795–804.
  * [26] D. Ionescu-Kruse,  On the short-wavelength stabilities of some geophysical flows,  _Phil. Trans. R. Soc. A_ , 376 (2018), 20170090, 21pp.
  * [27] M. Kluczek,  Exact and explicit internal equatorial water waves with underlying currents,  _J. Math. Fluid Mech._ , 19 (2017), 305–314.
  * [28] T. Lyons,  Particle trajectories in extreme Stokes waves over infinite depth,  _Disc. Contin. Dyn. Sys. Ser. A_ , 34 (2014), 3095–3107.
  * [29] T. Lyons,  The pressure distribution in extreme Stokes waves,  _Nonlin. Anal. Real World Appl._ , 31 (2016), 77–87.
  * [30] T. Lyons,  The pressure in a deep-water Stokes wave of greatest height,  _J. Math. Fluid Mech._ , 18 (2016), 209–218.
  * [31] T. Lyons,  The dynamic pressure in deep-water extreme Stokes waves,  _Phil. Trans. R. Soc. A_ , 376 (2018), 20170095, 13pp.
  * [32] C. I. Martin,  Dynamics of the thermocline in the equatorial region of the pacific ocean,  _J. Nonlin. Math. Phys._ , 22 (2015), 516–522.
  * [33] A. V. Matioc,  Exact geophysical waves in stratified fluids,  _Appl, Anal._ , 92 (2013), 2254–2261.
  * [34] W. J. M. Rankine,  On the exact form of waves near the surface of deep water,  _Philos. Trans. R. Soc. London A_ , 153 (1863), 127–138.
  * [35] A. Rodriguez-Sanjurjo,  Global diffeomorphism of the Lagrangian flow-map for Equatorially-trapped internal water waves,  _Nonlin. Anal.: Theor., Meth. Appl._ , 149 (2017), 156–164.
  * [36] A. Rodriguez-Sanjurjo and M. Kluczek,  Mean flow properties for equatorially trapped internal water wave–current interactions,  _Appl. Anal._ , 96 (2017), 2333–2345.
  * [37] K. W. S and M. M. J,  Oceaninc equatorial waves and the 1991–1993 El Niño,  _J. Climate_ , 8 (1995), 1757–1774.
  * [38] S. Sastre-Gomez,  Global diffeomorphism of the Lagrangian flow-map defining equatorially trapped water waves,  _Nonlin. Anal._ , 125 (2015), 725–731.
  * [39] G. G. Stokes,  Considerations relative to the greatest height of oscillatory irrotational waves which can be propagated without change of form,  _Mathematical and Physical Papers_ , 1 (1880), 225–228.
  * [40] J. F. Toland,  Stokes waves,  _Topol. Meth. Nonlin. Anal._ , 7 (1996), 1–48.
  * [41] G. K. Vallis, D. J. Parker and S. M. Tobias,  A simple system for moist convection: The rainy-Benard model,  _J. Fluid Mech._ , 862 (2019), 162–199.
  * [42] S. H. C. Wacogne,  _Dynamics of the Equatorial Undercurrent and Its Termination_ ,  PhD thesis, Massachusetts Institute of Technology, 1988.

Received August 2018; revised October 2018.