\section{Introduction}

The era of gravitational-wave (GW) astronomy began with the first direct
detection of the GW signal observed on September 14,
2015~\cite{GW170814-DETECTION} by the Advanced
LIGO~\cite{TheLIGOScientific:2014jea} and Virgo~\cite{Virgo-2015Acernese}
detectors. So far, in the proceeding observing runs, the LIGO Scientific and
Virgo Collaboration have fourteen confirmed GW
detections~\cite{GW150914-DETECTION,GW151226-DETECTION,GW170104-DETECTION,GW170608-DETECTION,GW170814-DETECTION,TheLIGOScientific:2016pea,O1-O2-catalog-paper,
GW170817-DETECTION, GW190425-DETECTION}, with data from several further
candidate compact binary merger candidates still being analyzed. In the next
few years, the Advanced LIGO and Virgo detectors will improve in sensitivity,
and additional GW detectors such as KAGRA \cite{Kagra-design}, and LIGO-India
\cite{LIGO-India} will join the network. The sensitivity improvement of the GW
detector network will lead to more detections of GWs from different types of
source~\cite{Davis:2020nyf}.

With improved detector perfromance, GW detectors are becoming more sensitive to
other disturbances culminating in the observation of a large number of
non-Gaussian noise transients known as {\it
glitches}~\cite{TheLIGOScientific:2016zmo}. They originate due to complex
instrumental and environmental effects within GW detectors. Often, they produce
high signal-to-noise ratio (SNR) triggers in GW searches, increasing false
alarm rates and reducing search sensitivity.
Data quality investigations~\cite{Nuttall_2015} that focus on the correlation
between instrument or environmental effects with GW data as well as veto
techniques have helped to eliminate some known glitch classes. However, a large
number of unknown noisy transients persist. Different searches explore
different types of glitch rejection
methods~\cite{Allen:2004gu,Nitz:2017lco,Canton:2013joa,Gayathri:2019omo} to
improve search sensitivity. Though important, these methods also increase the
computational burden of existing GW detection algorithms.

The volume of literature on the application of machine learning methods in GW
astrophsyics is rapdily growing. This includes techniques to improve signal
detection, parameter estimation of transient signals, noise removal techniques,
modeling GW signals as well as source population
inference~\cite{BAHAADINI2018172,PhysRevD.95.104059,Colgan:2019lyo,Cuoco:2020ogp,PhysRevD.101.042003,
George:2017pmj,Gebhard:2019ldz,Gabbard:2019rde,Cuoco:2020ogp}.

In this work, we propose a supervised machine learning method using {\it
Gaussian mixture models} (GMMs) for use in signal detection where we construct
two distinct models for noise glitches and astrophysical GW signals. Rather
than developing this approach using time-series data, we instead apply it after
the coherent Waveburst (cWB) algorithm has identified {\it triggers} --
interesting time instances which could be either potential GW signals or the
noisy transients. At present, the cWB search algorithm has multiple output
attributes and based on the location of the trigger attributes in this
multi-dimentional space, the trigger is classified as a GW event or noisy
event. We propose an alternative to the existing thresholding procedure on the
multi-dimentional attribute space by folding in the GMM naturally in the
detection problem under the likelihood ratio approach.

The paper is organized as follows : In Sec.~\ref{GMM}, we discuss the Gaussian
mixture modelling of the multi-modal data set. In Sec.~\ref{Sec:GMMLLR} we
discuss the use of the GMM in the construction of a log-likelihood based
detection statistic. In Sec.~\ref{burst} we assess the detection performance of the
proposed GMM based detection method with respect to the generic burst algorithm
in an all-sky short-duration GW burst set-up. We apply the algorithm for
the coincident events from the first observing run of advanced LIGO detectors.
Finally, in Sec.~\ref{conclusion}, we discuss our conclusions.

\section{Gaussian mixture model}
\label{GMM}
Gaussian mixture models (GMMs) are probabilistic models which use
uni-modal Gaussian distributions to represent a multi-modal data set. Under
the GMM approach, a given data set is modelled as a weighted sum of a
collection of Gaussians.

Let the data vector ${\bf x}$ be characterized by $d$ number of attributes. We
refer to these data as a $d$-dimensional data vector. The corresponding GMM of the data
consists of a superposition of $K$ Gaussian distributions
and is given by

\begin{align}
    p({\bf x}) = \sum_{j=1}^{K}w_j \mathcal{N}({\bf x}|\boldsymbol{\mu_j}, \Sigma_j),
\end{align}
where, $\mathcal{N} ( {\bf x}|\boldsymbol{\mu_j} , \Sigma_j )$ is a multinomial
Gaussian distribution with $d$-dimensional mean vector $\boldsymbol{\mu_j}$ and $d
\times d$ co-variance matrix $\Sigma_j$ and is written as
\begin{align}
   \mathcal{ N} ( {\bf x}|\boldsymbol{\mu_j} , \Sigma_j ) = \frac { \exp \left[ -\frac { 1
} { 2 } ({\bf x} - \boldsymbol{\mu_j} ) ^ { T } \Sigma_j^{ - 1 } ( {\bf x} - \boldsymbol{\mu_j})
\right] } { ( 2 \pi ) ^ { d / 2 } |\Sigma_j| ^ { 1 / 2 } }.
\end{align}
The parameter $w_j$ is the weight corresponding to each Gaussian component
normalised such that $\sum_j w_j=1$.

When the data is believed to contain two distinct populations but where each
population has complex but distinct structure within the attribute space, GMMs
can be used to model each population separately. These models can be then
incorporated into a likelihood-ratio test statistic which can be used to
identify events from either population. This is classed as an supervised
machine learning algorithm.

\section{Detection method using GMM}
\label{Sec:GMMLLR}
In the GW signal detection problem, signal events and noisy background
events are considered as two distinct populations. In this section, we describe
the use of GMMs to model these populations and the development of a detection
statistic based on this.

\subsection{Log-Likelihood statistic}
\label{sec:LLR}

Let us consider the data set with $n$ $d$-dimensional data points as ${\bf X} =
\{{\bf x}_1, {\bf x}_2, \dots{\bf x}_n\}$. We assume that each of the ${\bf
x}_i$'s are independent and we represent ${\bf X}$ as an $n \times d $ matrix.
The likelihood function can then be written as,
\begin{align}\label{maxL}
    p ( {\bf X}| \bm{\theta} ) = \prod _ { i = 1 } ^ { n } p ( {\bf x} _ { i }| \bm{\Theta}), \quad
\end{align}
where parameters $\Theta := w_j, \boldsymbol{\mu_j}, \Sigma_j, \{j=1,...,K \}$. Then the
corresponding total log-likelihood is the sum of $n$ individual log-likelihoods as given below,
\begin{equation}
    \ln \mathcal{L} =

    \sum _ { i = 1 } ^ { n } \ln ( p ( {\bf x} _ { i }| \bm{\Theta} ) )
    = \sum _ { i = 1 } ^ { n } \ln \left\{\sum _ { j = 1 } ^ { K } w_ { j }
\mathcal { N } \left( {\bf x} _ { i }| \boldsymbol{\mu_ j} , \Sigma _ { j }
\right)\right\} \, .\label{eq:like}
\end{equation}

\subsection{Maximum likelihood approach}
To estimate the model parameters $\Theta$  we can maximize  $\ln\mathcal{L}$
with respect to each of these parameters. For example,
\begin{align} \label{mu}
\frac { \partial \ln \mathcal{L} } { \partial   \boldsymbol{\mu_ j}}  = 0~~~ \Rightarrow~~~    \boldsymbol{\hat \mu_k} & = \frac { \sum _ { i = 1 } ^ { n } r _ { i k } {\bf x} _ { i }} { N _ { k } },
\end{align}
where
\begin{align}
r _ { i k } &= \frac { w_ { k } \mathcal { N } \left( {\bf x} _ { i }|
\boldsymbol{\mu_k} , \Sigma _ { k } \right) } { \sum _ { j = 1 } w_ { j }
\mathcal { N } \left( {\bf x} _ { i }| \boldsymbol{\mu_j} , \Sigma _ { j }
\right) }\hspace{0.2cm}\text{and}\hspace{0.2cm} N _ { k } = \sum _ { i = 1 } ^ { n } r _ { i k }.
\end{align}
Thus the maximum likelihood estimate of the mean $\boldsymbol{\mu_k}$ of the
$k^{\text{th}}$
Gaussian is the weighted mean of all the data points. All the coefficients are
implicit functions of ${\boldsymbol{\mu_k}}$ {\it via} the normal distributions.

Similarly, maximizing $\ln\mathcal{L}$ with respect to the co-variance matrix
$\Sigma_k$, of the $k^{th}$ Gaussian, we obtain
\begin{align}
    {\hat {\Sigma}} _ { k } = \frac { 1 } { N _ { k } } \sum _ { i = 1 } ^ { n } r _ { i k } \left( {\bf x} _ { i } - {\boldsymbol{\mu_k}} \right) \left( {\bf x} _ { i } - \boldsymbol{\mu_k} \right) ^ { T }.
\label{sigma}
\end{align}

Maximization of Eq.\ref{eq:like} over $w_j$ under the constrain that sum of the
weights add up to unity can be obtained through the application of Lagrange
multipliers. The details can be found in Appendix~\ref{ap1}. This gives the maximum likelihood estimate for weights as
\begin{align}
\hat{w}_ { k } &=  \frac{\sum _ { i = 1 } ^ { n } r_{ik}} {n}  = \frac { N _ {
k } } { n }.\label{maxW}
\end{align}

It is clear from above calculation that it is difficult to analytically
estimate the parameters of the mixture model as all the estimates given in
Eq.~\ref{maxW} are implicit functions of themselves.

{The \it{Expectation maximization} (EM)
technique}~\cite{doi:10.1111/j.2517-6161.1977.tb01600.x} is an iterative
algorithm and provides us with a numerical solution to our maximum likelihood  problem.
The EM algorithm has two steps, namely the estimation step and the maximization step. The first
step involves using trial values for the parameters, then, using these values
an iterative step using Eqs.~(\ref{mu}, \ref{sigma}, \ref{maxW}) gives
estimates of the new values of the parameters. Thus, iteratively using
Eqs.~(\ref{mu}, \ref{sigma}, \ref{maxW}), the EM algorithm convergences on the
maximum-likelihood parameters $\hat{\Theta}$.

A drawback of the EM algorithm is that it cannot predict the optimal number of
Gaussians required to describe the underlying structure of the data. Large
numbers of Gaussians can lead to overfitting of the data. The {\it{ Bayesian
information criterion}} (BIC) includes a penalty term which compensates this
effect and is used in the model selection.

The BIC as defined in terms of the maximum value of the likelihood function
$\hat{\mathcal{L}}$ and is given by

\begin{align}
    {\rm BIC} &= K\ln(n) - 2\ln(\hat{\mathcal{L}}).\label{eq:BIC}
\end{align}
where the first term in Eq.~\ref{eq:BIC} is the desired penalty term. The lowest
BIC score provides the optimum number of Gaussians $\hat{K}$ for a given data set.

\subsection{GMM based detection statistic}

Once all the model parameters $\Theta$ are optimally chosen and the optimum number of Gaussians are fixed following the minimum BIC criterion detailed in the previous subsection, we
write the maximum log-likelihood statistic as $W = \ln(\hat{\mathcal{L}} )|_{\hat{K}}$.

Since our data consists of two distinct classes, signals (s) and noisy background glitches (g), we can calculate a detection statistic, $T$, for each trigger such that
\begin{align}
\label{T}
T = W_{\text{s}}- W_{\text{g}}\,.
\end{align}

GMM based models each for the signal and noise are required to determine $W_{\text{s}}$
and $W_{\text{g}}$ respectively.
In line with tuning and training procedures for cWB and other transient searches, the characterisation and optimisation of the GMM performance is done {\it a prior} on training data, prior to the search for gravitational wave candidates. To calculate $W_{\text{g}}$, we first construct a GMM model using the noise background data.
The noise background data are divided into a training data set, which is used to construct the GMM model, and a validation data set.
Similarly, for $W_{\text{s}}$, we construct a GMM model using simulated signals which are also divided into training and validation sets.
We compute the test statistic $T$ for the validation noise set and validation simulation trigger set and we assess the performance of the GMM models by comparing the detection probability against the false alarm probability. Note that the validation data set allows us to check that the GMM model is overfitting the signal parameter space since it is vital for a generic transient search to be sensitive to a wide range of signal morphologies.

\section{Short duration GW burst search with GMM}

\label{burst}
A large variety of GW transients fall into the category of short duration
bursts e.g., GWs from supernovae, merger of binary black holes etc. In fact,
the first observed GW signal from a binary black hole merger (GW150914) was
indeed first detected by a generic burst search algorithm coherent Waveburst
(cWB)~\cite{GW150914-DETECTION}. Typically, the cWB burst search method finds
interesting and potential GW events at an initial analysis stage which we refer
as {\it triggers}.  To obtain an estimate of the noise background and assign
each trigger a statistical significance, the cWB analysis is performed on data
that is time-shifted so that it is unphysical for a gravitational wave signal
be detected in coincidence between detectors. The noise background rate is
determined by the number of triggers observed in the time-shifted data.
Triggers observed in data in the absence of an unphysical time shift are
considered {\it event candidates}. These event candidates may be GW signals
though further analysis is often required before any declaration of signal
detection can be made.

To reduce the impact of noise background triggers, thresholds for various
trigger attributes such as its estimated signal strength or duration are
applied. To optimise the chance of detecting a GW signal, these thresholds are
tuned {\it a priori}, using triggers from the time-shifted noise background and simulated
signals. After applying appropriate thresholds on the trigger attributes
pertaining to the search, if a given trigger emerges with high significance
then it is considered a GW event. Such a thresholding procedure, though ad-hoc,
is the crucial ingredient of the signal detection algorithm which helps in
reducing the unwanted noisy features and select the most significant events
(those with very low associated false alarm rate (FAR)) in the data set.

Here, we introduce GMMs to address this multi-dimensional attribute
thresholding problem where a GMM is constructed for the signal set (and noise
set) defining the characteristic features in the multi-dimentional attribute
space for signal (or noise). Casting this in the log-likelihood based detection
problem as outlined in Sec.~\ref{Sec:GMMLLR} gives a natural map from this
multi-dimentional approach to the scalar log-likelihood ratio. Thus, by
applying a single threshold based on the test statistic across the entire
multi-dimensional attribute space, we can address the ad-hoc individual
attribute thresholding problem in a more systematic way. We henceforth refer to
this signal detection approach as cWB plus GMM.

\subsection{Coherent Waveburst algorithm}
\label{cWB}
For our data sets we consider burst-search triggers from the coherent waveburst
(cWB) algorithm -- a generic multi-detector, all sky burst algorithm used in GW
searches that is sensitive to short duration
signals~\cite{Klimenko:2008fu,klimenko05:_const,Klimenko:2015ypf,
necula12:_trans_wilson_daubec}. This algorithm projects the multi-detector data
into the wavelet (time-frequency map) domain using the
Wilson-Daubechiers-Meyer transformation and identifies a collection of coherent
time-frequency-scale pixels with excess power and clusters them based on the
time-frequency information. Each burst trigger is represented by a coherent
cluster and associated with a set of attributes. The attribute set contains those which characterise
the signal as well as veto attributes used to help distinguish between signal and noise triggers. The complete set include estimated strain of the trigger $h_s$, the central frequency
$f_0$, the duration of the trigger $\tau$, the network coherent signal-to-noise
ratio $\eta_c$~, the network correlation coefficient $c_c$ which measures
the correlation between the detectors,
quality factor of the event $Q_{\text{veto}}$, the residual noise energy measure
$\zeta^2$, the ratio between the reconstructed energy and the total energy
$N_{\text{norm}}$, the energy dis-balance of the event between the detectors
$N_{ED}$, and $L_{\text{veto}}$ which measures the localization of event in the
time-frequency map.

The cWB algorithm generates a trigger list that is associated to the
time-frequency clusters with high network correlation $c_c$ as well as cWB SNR
$\eta_c$. Based on the type of the GW signal, the actual threshold values
applied to these two parameters may vary. Following that, the cWB algorithm
uses the additional thresholds on trigger attributes to veto out noisy
triggers. All cWB thresholds are determined by tuning and testing on noise background and simulated signals prior before being used to generate a list of gravitational wave candidates.

\subsection{Data set}
\label{sec:IVB}

To demonstrate our method, we consider cWB triggers obtained from simulated
astrophysical short-duration burst signals added to data from the first advanced
detector observing run O1~\cite{Abbott:2016ezn}. Here, cWB triggers are those
time-frequency clusters with $\eta_c >8.45$ and $c_c > 0.5$.  The 19
types of short duration signals simulated in the O1 data are {\it Gaussian pulses} (GP) characterized by duration $\tau$, {\it sine-Gaussian
wavelets} (SGW) - sinusoids within a Gaussian envelope and characterized by the
frequency$f_0$  and a quality factor $Q$, and {\it White noise bursts} (WNB) -
bursts with a Gaussian envelope described by $f_{\text{low}}$, frequency
bandwidth $\Delta f$, and duration $\tau$. The simulated signal parameter
values used our study is listed in Table~\ref{Tab:BURST}. The signals are
uniformly distributed over the sky and with a range of selected discrete
strain amplitudes from $0.5 \times10^{-21}$-$10^{-20}$. The initial phase is
distributed uniformly over the range $[0,2\pi]$ and the time of arrival of the
signal which depends on the sky location in relation to the detector positions is also distributed uniformly.

The final GW event list is obtained with the standard-cWB after applying
thresholds on the cWB attributes and ranking them based on level of
significance (ascending value of inverse false alarm rate (IFAR)). For more
details on the search see~\cite{Abbott:2016ezn}.

\subsection{Attribute choice}
\label{sec:opt-para}

Each trigger has eleven attributes as defined in Sec.~\ref{sec:IVB}. Here, we
use the exhaustive attribute set $\{c_{c0},c_{c2},\eta_c,f_0, h_s,
N_{\text{norm}},\zeta^2 ,N_{\text{ED}},Q_{\text{veto}},L_{\text{veto}}\}$ to
develop a GMM in the multi-dimentional feature attribute space and use the
minimum BIC method to obtain the optimal attribute set which we use to
construct the GMM for both noise as well as signal data triggers.

For a collection of $n \sim 10^5$ signal triggers from the data set, we
consider the maximum number of Gaussian components to be $100$. We carry out an
exhaustive study to build the signal model. We consider various
combinations of trigger attributes starting from a set of two attributes all
the way to a set of eleven attributes with a varying number of Gaussian
components. We compute the lowest BIC score for a given number of attributes.
For example, there are 55 combinations of pairs of attributes selected from a
set of 11, 165 sets of 3 attributes selected from 11, etc. We construct the BIC
for all possible combinations and choose the minimum value amongst them for a
given attribute count.

This is done for each attribute count ranging from 2 to 11.

\begin{figure}[ht]

    \centering

    \begin{minipage}[b]{.45\textwidth}
        \includegraphics[scale=0.35]{minimum_bic_sim.png}
    \end{minipage}
    \caption{Plot of BIC value {\it{vs}} the number attributes for a given signal data
set. The specific choice of attributes for each set is defined in the
legend.
}
    \label{fig:min_bic}
\end{figure}

\begin{figure*}[!ht]
\includegraphics[scale=0.4]{P_dist_model_all_model_paper_zoom_T5_nofreq_july15_steps.png}
\caption{The $T$ value distribution for O1 noise and three  different  signal
morphology events with GMM models. Plot contains four curves,  red curve for
noise, blue curve for GP, green curve for SWG and orange curve for WNB. The
inset plot shows the T distribution around zero.}
\label{fig:F2}
\end{figure*}

\begin{figure*}[!ht]
    \hspace{-0.5 cm}
    \centering
    \begin{minipage}[b]{.45\textwidth}
        \includegraphics[scale=0.25]{1_withshift_with_cwb_july10.png}
    \end{minipage}
    \qquad
    \begin{minipage}[b]{.45\textwidth}
        \includegraphics[scale=0.25]{2_withshift_with_cwb_july10.png}
    \end{minipage}
    \hfill
        \begin{minipage}[b]{.45\textwidth}
        \includegraphics[scale=0.25]{3_withshift_with_cwb_july10.png}
    \end{minipage}
    \qquad
    \begin{minipage}[b]{.45\textwidth}
        \includegraphics[scale=0.25]{4_withshift_with_cwb_july10.png}
    \end{minipage}
    \caption{ROC curves: false alarm probability {\it vs} detection probability using the cWB plus GMM algorithm (dashed lines)
compared to the standard-cWB pipeline (solid lines) for the simulation waveforms detailed in Table.\ref{Tab:BURST}. We present the ROC curves in four panels by dividing the simulation sets in different types. For the cWB plus GMM
analysis, the detection probability is close to 0.9 or more for almost all
simulated signals, even at false alarm probabilities of $10^{-6}$.}
\label{fig:ROC}
\end{figure*}

Figure~\ref{fig:min_bic} shows a plot of the minimised BIC value as a function
of the attribute count for the signal model. The combination of attributes
shown in the legend corresponds to the minimum BIC combination for that
attribute count. We note that the minima of this minimum BIC curve corresponds
to the attribute set [$\zeta^2$, $N_{\text{norm}}$, $h_s$, $\tau$, $c_{c0}$,
$c_{c2}$,\footnote{Both $c_{c0}$ and $c_{c2}$ are network correlation coefficients and capture signal correlation between the detectors.} $Q_{\text{veto}}$]-- a seven attribute set which captures the signal triggers. In addition, we note that each set of $p$ optimal attributes contains within it the set of $(p-1)$ optimal attribute set.
We use the attribute set corresponding to the global BIC minima to construct both the signal and noise trigger
GMMs.

\subsection{Training models}

As described in Sec.~\ref{sec:IVB}, the simulated signal injections span a
variety of signal classes that model short bursts with a wide range of
frequencies.
We aim to develop a GMM for these signals that is robust against the variation
in frequency and duration of the burst. To do so, we develop a minimal model
which broadly captures the burst-like signals with a frequency range of
$70-360$ Hz and different short durations.
Thus, we choose a waveform subset from  the simulation set given in table. \ref{Tab:BURST} which cover the broad class of burst-like signals for model building and test the robustness of this model against the remaining waveform set.  With
the above rationale, our choice is a set of five burst waveforms; SGW with
$f_0=70$Hz \& $Q=3$, SGW with $f_0=235$Hz \& $Q=9$, SGW with $f_0=361$Hz \&
$Q=9$, SGW:$f_0=70$Hz \& $Q=100$ and GP $\tau=1$s as tabulated in the
Table~\ref{Tab:BURST}.

We choose 50\% of the five selected types of burst signal triggers for
training and building the GMM. The remaining signal triggers (comprising all
waveform types) are used for validation, to check that the signal model is robust against varying signal morphologies.
We consider 80\% of noise triggers for
training and building the corresponding noise GMM.
The rest of the noise triggers are used for validation.  This leads to the number
of noise triggers for training being $\sim 10^5$ and signal triggers for training
as $\sim 8 \times 10^4$.

\begin{table}[htb]
\begin{tabular}{c c c c c c}
\hline
\hline
\multicolumn{5}{c}{Sine-Gaussian Burst (SGW)}\\
\hline
No. &$f_0$ (Hz) & $Q$ & - & Training & Validation\\
\hline
1 & 70 & 3 & - & Y & Y\\
2 &70 & 9 & - & N & Y\\
3 &70 & 100 & - & Y & Y\\
4 &100 & 3 & - & N & Y\\
5 &153 & 9 & - & N & Y\\
6 &235 & 3 & - & N & Y\\
7 &235 & 9 & - & Y & Y\\
8 &235 & 100 & - & N & Y\\
9 &361 & 9 & - & Y & Y\\
10 &554 & 9 & - & N & Y\\
11 &849 & 3 & - & N & Y\\
12 &849 & 9 & - & N & Y\\
13 & 849 & 100 & - & N & Y\\
\hline

\multicolumn{5}{c}{White-Noise Burst (WNB)}\\
\hline
 & $f_{low}$ (Hz) & $\Delta f$ (Hz) & $\tau$ (s)& Training & Validation\\
\hline
14 & 100 & 100 & 0.1 & N & Y\\
15 & 300 & 100 & 0.1 & N & Y\\
\hline

\multicolumn{5}{c}{Gaussian Pulse (GP)}\\
\hline
 & - & -& $\tau$ (s) & Training & Validation\\
\hline
16 & - & -& 0.1 & N & Y\\
17 & - & -& 1 & Y & Y\\

18 & - & -& 2.5 & N & Y\\
19 & - & -& 4 & N & Y\\
\hline
\hline
\end{tabular}
\caption{List of generic burst waveforms and their characteristic parameters using in training and validation the GMM signal model.}\label{Tab:BURST}
\end{table}

\subsection{Search sensitivity improvement}
\label{Sec:SensImprove}

In this section, we assess the improvement in search sensitivity by comparing
signal recovery with the standard-cWB algorithm and the cWB plus GMM based
detection algorithm.

After constructing the GMMs for signal and for noise triggers,
Fig.~\ref{fig:F2} shows the distribution of the detection statistic $T$ for
noise and signal testing data computed using Eq.~\ref{T}. The plot shows four
curves, noise (red), GP signals (blue), SWG signals (green), and WNB signals
(orange). We observe that the majority of noise triggers have negative $T$
value and similarly most signal triggers have positive $T$ value. There is a
corresponding clear separation between the signal and the noise triggers at the
$T=0$ boundary. Very small overlap exists between the T distributions for
signal and noise triggers. This clearly indicates that the GMM  based proposed
likelihood ratio statistic has the potential to discriminate between the two
classes of triggers.

To compare the performance between the standard-cWB and the cWB plus GMM based
detection approach, we construct Receiver Operating Characteristic (ROC)
curves. ROC curves relate the detection probability and the false alarm
probability as the threshold on the test statistic is changed. For the
standard-cWB analysis, this means recording the fraction of simulated signals
and noise triggers for varying threshold values on $\eta_c$, while for cWB plus
GMM, it is the threshold on $T$ that is varied. We use $\sim 3 \times 10^4$
noise triggers and $10^5$ signal triggers to create the ROC curves shown in
Fig.~\ref{fig:ROC}.

We note that cWB plus GMM based detection algorithm significantly enhances the
signal detection probability when compared to the standard-cWB detection alone.
The detection probability using the cWB plus GMM varies between 0.9 and 1 for
most of the waveform models. The performance for SWG waveform with
$f_0=849$Hz and $Q=100$.  (see panel c in Fig.~\ref{fig:ROC}) appears to have dropped a
little compared to other waveforms. This is primarily because the cWB plus GMM
approach is sensitive to the parameter distributions used within the training
data and we have chosen to use burst waveforms with frequency values up to 350
Hz and short durations. For high frequency signals, although we do not get as spectacular
performance as the other cases, we do achieve significantly improved performance
compared to the standard-cWB settings.

\subsection{Application of cWB plus GMM to the O1 data}
\label{Sec:AppToO1}

We apply the cWB plus GMM algorithm to the coincident data from O1. The GMM
signal model is already trained using the all-sky short duration cWB triggers obtained
from simulations as detailed in Sec.\ref{Sec:SensImprove}. The GMM noise model is also
trained on the time-shifted cWB background triggers from the O1 all-sky short
duration search. We compare the results after applying the cWB plus GMM
algorithm with the standard-cWB results obtained by the all sky short duration
burst search reported in \cite{Abbott:2016ezn}.

For the O1 all-sky short duration search, a total of 10 cWB triggers are
observed when no time shift is applied to the data.  For the standard-cWB
algorithm, additional thresholds are applied on the cWB trigger attributes to
further reject the noise triggers. The additional threshold requirements used
in the standard-cWB used for the O1 all-sky short duration search are $c_{c0} >
0.7$, $48$Hz$<f_0<998$Hz, $Q_{\text{veto}}>0.3$ and $\zeta^2 <0.5$. A total of
6 coincident events survive these threshold criteria~\cite{Abbott:2016ezn}. We
draw these events as red hollow circles in Fig.~\ref{fig:zero}.

Instead of applying various thresholds on individual trigger attributes, we
consider the cWB plus GMM based detection approach on the triggers. We compute
the $T$ statistic for each cWB trigger and select events with $T>2.5$
(corresponding to an IFAR of 15 years) as candidate events denoted by
large black hollow circles in Fig.~\ref{fig:zero}. For the IFAR estimation, we consider the entire background trigger set.

Figure~\ref{fig:zero} shows $\eta_c$ vs $T$, with color corresponding to
$c_{c0}$ for all cWB triggers ($\eta_c>8.45$ and $c_{c0}>0.5$). The inset plot
shows $\eta_c$ vs $T$ for the events clustered below $\eta_c=8.9$. We show all
10 cWB triggers as filled circles with the colour indicating the value of
$c_{c0}$. We notice that 5 events are common in both the searches. In addition,
the cWB plus GMM identifies three events which were rejected by the
standard-cWB analysis and one event from standard-cWB is rejected by cWB plus
GMM.

Both standard-cWB and cWB plus GMM observe GW150914 to be the most significant
event, which is indicated by a star in Fig.~\ref{fig:zero}. For cWB plus GMM,
GW150914 is observed with $T\sim41$. With no background events having that $T$
value or greater, the cWB plus GMM analysis detects GW150914 with an IFAR of
greater than $1000$ years. We note that the standard-cWB analysis for O1
detects GW150914 with an IFAR of $\sim 350$ years, even after the analysis of
O1 data was split into 3 classes to isolate chirp-like signals into a single
class to separate them from spurious noise transients~\cite{Abbott:2016ezn}. The remaining four common events also show higher significance in cWB plus GMM
compared to the standard-cWB. The two additional events observed by the cWB
plus GMM have network correlation just below 0.7 and hence went undetected by
the standard cWB.

\begin{figure}[t]
\includegraphics[scale=0.4]{paper_eta_c_T_vs_cc_paper_nofreq_loglog_v3.png}
\caption{The network coherent SNR $\eta_c$ vs the GMM detection statistic $T$,
with color corresponding to $c_{c0}$ for all 10 cWB triggers.
Large black hollow circles correspond to detected cWB with GMM events (with
$T>2.5$ equivelent to IFAR$>15$years) and smaller red hollow circles correspond to the
standard-cWB events. The black star corresponds to the GW150914 event}
\label{fig:zero}
\end{figure}

\section{Conclusion}
\label{conclusion}

Detection of transient GWs in non-stationary and non-Gaussian noise poses a
massive challenge in the advanced laser interferometric gravitational wave
detectors. A number of approaches are used to combat spurious noisy triggers
and to allow the detection of astrophysical GW transients with high
significance.  The standard-cWB GW burst search {\it first} makes a list of
triggers -- interesting time instances which could be either potential GW
signals or background noise transients -- and {\it then} applies a list of
thresholds on the attributes to reject the noisy triggers.  Although this
thresholding approach is guided by extensive simulations and hence is largely
successful, it is still ad-hoc in nature and relies on analysts intuition to
optimise the thresholding process. In this work, we propose a supervised
machine learning method using GMMs trained on cWB triggers.  This training is
applied using the O1 all-sky search triggers from the cWB search for short
duration GW bursts to develop two distinct models for noise triggers and
astrophysical GW burst triggers. We use these models to construct a
log-likelihood based test statistic. We demonstrate that this approach gives
improved performance as compared to the standard-cWB approach by improved
signal detection probability at any FAR. We also obtain a significantly
improved detection significance for GW150914 (the first ever GW merger event)
with the cWB plus GMM detection approach. With this example, we clearly
demonstrate that our more systematic GMM based signal detection approach can
improve the detection performance as compared to the thresholding approach in
the multi-dimensional attribute space.

It is worth noting that, like most GW search algorithms, the GMM approach is
sensitive to simulation parameters used to train the signal model. Here, for
comparison with standard-cWB, we have used the population of simulated signals
that were generated for the O1 all-sky short duration burst search. The
simulated data set was designed to provide an estimate of the search
sensitivity to an ad-hoc population of signals and was not drawn from an
astrophysical population. For example, the simulated data did not uniformly
span the frequency band of the detectors. If we had included frequency as one
of the parameters for the GMM, then the resulting test statistic would favour
signals with frequencies corresponding to the simulated population which is not
a desirable feature for unmodelled burst searches.

By applying the GMM to cWB trigger attributes, we are relying on cWB to
efficiently capture all characteristics of an interesting event in these
attributes. The cWB algorithm has undergone almost two decades of development
and much effort has been put into the development and characterisation of the
algorithm. Nonetheless, it may be possible to find more efficient ways to
encode GW data for more effective GMM signal and noise model construction. For
example, neural networks can be used to map GW strain data into a reduced
number of parameters which are then used as inputs into the corresponding GMM.

In the near term, we plan to test this GMM approach using data from the second
advanced detector observing run which contains a large variety of noisy
triggers. The method is general enough and is not limited to short duration
bursts triggers. We plan to extend this approach to GW searches for
intermediate mass binary black hole mergers as such signals can be very
difficult to distinguish from spurious noise transients, especially with
increasing black hole mass.

\section{Acknowledgements}
The authors thank F. Salemi for useful comments. VG acknowledges Inspire division, DST, Government
of India for the fellowship support. VG acknowledges support from the University of Florida. DL and RP thank Newton-Bhabha fund for the travel
support to visit the University of Glasgow. AP acknowledges support from the SERB
Matrics grant MTR/2019/001096 as well as the IIT-Bombay SEED grant for the
travel funds to host the visit of ISH and CM to IIT Bombay. CM and ISH are
supported by the Science and Technology Research Council (grant No. ST/
L000946/1) and the European Cooperation in Science and Technology (COST) action
CA17137.


\section{Introduction}\label{sec:introduction}
Localization phenomena play a predominant role in Earth Sciences as many interesting geological features fall within this category, including faults, folds, boudinage, landslides, and mineralization, just to name a few. Among all localized geological features, spatially periodic patterns are increasingly gaining popularity due to their impact in a wide range of applications, such as mining~\cite{Iophis2007}, particularly now when larger-scale and deeper mines are pursued. Localization features are particularly relevant when they affect permeability, which plays a critical role in various fields including mineral exploration~\cite{Hayward2018, Kelka2017} and unconventional resources exploration~\cite{Regenauer-Lieb2016} as well as their exploitation.

One mechanism responsible for spatially periodic patterns affecting porosity and permeability is the formation of compaction bands. These bands are narrow flat zones of deformation perpendicular to the maximum compressive principal stress~\cite{Das2013}. For example, a succession of compacted zones of lower permeability in a non-compacted higher permeability background is the most intuitive way to imagine compaction bands~\cite{Holcomb2003}, with critical implications on the fluid flow in terms of creation of impermeable barriers as well as compartmentalized reservoirs and aquifers. These compaction bands therefore have significant impact on the fluid production or geological storage (${\text{CO}_2}$, nuclear waste). Several models exist to describe the mechanism on its in situ periodic occurrence, e.g.~by Cecinato~\cite{Cecinato2014}, however, there also exist other physical processes that lead to periodic bands under compression with increased permeability, an example from the melt segregation field known as decompaction bands~\cite{Rabinowicz2004}.

Attempts to theoretically describe deformation bands (principally of shear-type) start early in the theory of plasticity~\cite{Hill1950}. For geomaterials, this analysis started with the work of Rudnicki and Rice~\cite{Rudnicki1975}, who extensively study the conditions for localized deformation in low porosity rocks using bifurcation analysis. Although that theory sought to explain the appearance of shear bands in a specific type of rocks, a subsequent work by Olsson~\cite{Olsson1999} confirmed that the origin of compaction bands could also be explained within this framework. Later, a re-examination of the original theory was developed by Issen and Rudnicki~\cite{Issen2001}, where the triggering conditions of both shear and compaction bands could be deduced for a broader range of materials. Thus, they proposed the introduction of a cap-type yield function to take into account compactive yield in high porosity rocks. Finally, in the particular case of cemented soils, studies on the onset of compaction bands were developed by Arroyo et al.~\cite{Arroyo2005} using a bonded soil model and bifurcation conditions established by Rudnicki and Rice~\cite{Rudnicki1975}.

From a numerical perspective, little has been done on the simulation of compaction bands compared with the numerous studies on the formation of shear banding~\cite{Oka2011}. Within the most relevant works can be named the computational modeling carried out by Borja~\cite{borja2004}, where the onset conditions for deformation bands were analyzed using the single hardening constitutive model proposed by Kim and Lade~\cite{Kim1988}. In addition to this, Oka et al.~\cite{Oka2011} developed finite element simulations in diatomaceous mudstone using an elasto-viscoplastic model and compared them with triaxial tests results. Although their models predicted compaction bands for higher confining pressures, they did not tackle the challenging problems of identifying the onset conditions of the phenomenon and its periodicity for a broader range of confinement stresses. Other experimental studies were also developed for sandstone specimens~\cite{Fortin2006} and calcarenite samples~\cite{Baxevanis2006}.

All previous studies only focused on the mechanical problem, without considering multi-physical scenarios that could allow capturing other contributions to the deformation process. To overcome this issue, recently, Veveakis and Regenauer-Lieb~\cite{Veveakis2015} developed a wave mechanics approach which showed that regularly spaced localization bands of hydro-mechanical nature can appear in rocks under compaction. The wave description is radically different as this type of localized deformation stems from a generalization of Terzaghi's linear theory of consolidation to materials with non-linear visco-plastic rheology and any arbitrary type of internal mass transfer mechanism. Through these genearlizations, the authors derived the following governing equation for nonlinear consolidation, admitting material instabilities presented in the effective stress $\sigma’$ (for elasto-viscoplastic materials subject to hydro-mechanical processes):

\begin{equation}\label{eq:cnoidal_orig}
\frac{\partial^2\sigma'}{\partial{z^2}} - \lambda\sigma'^m =0
\end{equation}
In this dimensionless equation $\lambda$ represents a ratio between the particular diffusive processes of the problem (i.e., the mechanical deformation of the matrix and the internal mass exchange) and $m$ is a material-dependent pressure exponent. Under certain conditions, the solution of~\eqref{eq:cnoidal_orig} presents numerical instabilities due to the loading rate being faster than the mass diffusion rate; thus mass variations in the specimen cannot be equilibrated, producing zones of stress concentrations that represent the compaction bands. This phenomenon produces periodic volumetric failure patterns denoted as cnoidal waves~\cite{Regenauer-Lieb2016,Veveakis2015,Veveakis2014}, due to the analogy with its counterpart in fluids dynamics. “Cnoidal wave” is the term used for the solution of the Korteweg-de Vries (KdV) equation~\cite{Korteweg1895} expressed regarding the square Jacobi $cn$ elliptic function. The KdV equation describes a traveling wave in shallow water surfaces, first observed by Russell~\cite{Russell1844}, which has been extensively used for many physical problems related to wave mechanics. For the case of solids, this cnoidal wave approach allows us to conclude that periodic instabilities are produced as a volumetric response during failure and that the latter is controlled by deformation rate, and not by critical stress or hardening reaching as stated by classical theories~\cite{Veveakis2015}.

While the original theory~\cite{Veveakis2015} captures well the essence of hydro-mechanical instabilities, its simplifying nature leads to unbounded stress values, which are not realistically possible. A complementary study~\cite{Alevizos2017} considered, for instance, the effects of chemical reactions, which introduce a regularization term in the equation and allow the effective stress to remain capped. Without entering the debate about which physical processes could occur on the back of stress peaks, we find if convenient to consider a nonlinear regularization term $N_r(\sigma')$ for mathematical reasons. This regularization allows us to deal with bounded continuous solutions. As such, we consider a generalization of equation~\eqref{eq:cnoidal_orig} of the form
\begin{equation}\label{eq:cnoidal_chem}
\frac{\partial^2\sigma'}{\partial{z^2}} - \lambda\sigma'^m + N_r(\sigma') = 0,
\end{equation}
and for this study, we use the existing formulation of~\cite{Alevizos2017}. We succinctly recall it in Appendix~\ref{AppxA} for completeness.

In this work, we extend an adaptive stabilized finite element method (FEM) to overcome the numerical problems by which stress singularities hinder the effectiveness and robustness of numerical schemes. We showcase the power of such approach by solving this equation consistently. Next, we implement an arc-length continuation algorithm to perform a numerical bifurcation analysis to better understand the parameter sensitivity of the model. This analysis allows us to derive important conclusions regarding the stability regimes of the system. We introduce the numerical approach used for the resolution of~\eqref{eq:cnoidal_chem} in Section~\ref{sec:num} and all the numerical outcomes in Section~\ref{sec:numexp_1d}, leading to the discussion in Section~\ref{sec:discussion}.

\section{Numerical approximation of the cnoidal problem}\label{sec:num}
The cnoidal wave approach in solids seeks to explain the formation of specific localized deformation bands more readily than alternative explanations provided by classical theories. While the cnoidal approach offers a new perspective to the localization phenomenon, there are some points related to the solution of its governing equation that need to be addressed before a detailed study is possible. Equation~\eqref{eq:cnoidal_chem} has known analytical solutions only for the integer values of $m=1, 2, 3$. However, solutions for higher or non-integer values of $m$ need numerical treatment, and, to date, there has been no successful attempt to numerically solve this equation satisfactorily. The lack of a robust numerical solution for this equation is related to the complexity of the treatment of this class of nonlinear problems. We seek to overcome this issue by developing a consistent numerical solution for~\eqref{eq:cnoidal_chem}. Therefore, we use the new adaptive stabilized finite element method based on residual minimization, developed by Calo et al.~\cite{calo2019}, and extend its application to nonlinear problems. We opt for using this formulation instead of rather than alternative FEM approximations because of the stability properties this method enjoys, and its built-in adaptive mesh refinements, a crucial feature for localization of instabilities.

\subsection{Weak variational formulation}
In an abstract setting, we consider a well-posed variational formulation for a general nonlinear problem. For an open set $\Omega\neq0$ with boundary $\partial \Omega$, and Hilbert spaces $U$ (trial) and $V$ (test), let $N$ be a differentiable nonlinear map with Fr\'echet derivative $DN(u)$ at ${u \in \Omega}$. We associate the nonlinear map ${n:U\times V \rightarrow \mathbb{R}, \, n(u;v):=\langle N(u),v \rangle}$, where ${\langle \boldsymbol{\cdot}, \boldsymbol{\cdot} \rangle}$ represents the duality pairing in $V$. Let ${n'(u;z,v)}$ abbreviate the derivative $DN(u)$, around a known value $u$, and in the direction of an increment $z$:
\begin{equation}\label{eq:der-def}
 n'(u;z,v):=\langle DN(u;z),v\rangle = \frac{d}{d\epsilon}n(u+\epsilon z; v)\big|_{\epsilon=0} \quad \text{for } u \in \Omega, z \in U, v \in V.
\end{equation}
We finally set $\ell(\boldsymbol{\cdot}): V \rightarrow \mathbb{R}$ as a continuous linear form. Hence, the weak formulation for a nonlinear problem reads:
\begin{equation}\label{eq:abs_vf}
\left\{
\begin{array}{l}
\text{Find } u \in U, \text{ such that:} \smallskip \\
\begin{array}{l}
n(u;v)= \ell(v), \quad \forall v \in V,
\end{array}
\end{array}
\right.
\end{equation}

We now state the variational formulation for the particular case of the cnoidal equation~\eqref{eq:cnoidal_chem}. Following~\cite{Veveakis2015}, we define the boundary conditions as $\sigma' = 1 $ on $\partial \Omega$. In order to derive the continuous formulation of~\eqref{eq:cnoidal_chem}, we split ${\sigma'=u+1}$, and the steady state equation with homogeneous Dirichlet boundary conditions then reads:
\begin{equation}\label{eq:numsol_def_equation}
\left\{\begin{array}{l}
\text{Find } u=\sigma' -1 \text{ such that:} \smallskip \\
\begin{array}{rcl}
\Delta u-\mathcal{F}(u)&=&0 \quad \text{in} \quad \Omega,\\
u&=&0 \quad \text{on} \quad \partial \Omega,
\end{array}
\end{array}\right.
\end{equation}
with ${\mathcal{F}(u)=\lambda\, (1+u)^m -\mu \exp(\beta u)}$. In~\eqref{eq:numsol_def_equation}, ${\mu \exp(\beta u)}$ is equivalent to the regularization term $N(\sigma')$ from~\eqref{eq:cnoidal_chem}, which avoids the unbounded stress growth of the exponent ${m>1}$ capping the value of $u$ to finite values and making $\mathcal{F}(u)$ to remain positive.

\sloppy We consider the standard notation for Hilbert spaces: ${L^2(\Omega)=\{w: \Omega\in \mathbb{R}^d  \rightarrow \mathbb{R} : \int_{\Omega}|w|^2<+\infty\}}$,  ${H^1(\Omega)=\{w\in L^2(\Omega): \nabla w \in [L^2(\Omega)]^d\}}$ and ${H^1_0(\Omega)=\{w\in H^1(\Omega): w=0 \; \text{on} \; \partial\Omega\}}$. Multiplying~\eqref{eq:numsol_def_equation} by a test function $v$ and integrating by parts, we obtain the following weak variational formulation:
\begin{equation}\label{eq:weak_nonlin}
\left\{\begin{array}{l}
\text{Find } u \in H^1_0(\Omega), \text{ such that:} \smallskip \\
\begin{array}{l}
n(u;v)= \ell(v), \quad \forall v \in H^1_0(\Omega),
\end{array}
\end{array}\right.
\end{equation}
with $n(u;v)=(\nabla u, \nabla v)_{0, \Omega} + (\mathcal{F}(u), v)_{0, \Omega}$ and $\ell(v)=(f,v)_{0, \Omega}$, where $(\boldsymbol{\cdot}, \boldsymbol{\cdot})_{0, \Omega}$ represents the $L^2$ scalar product in $\Omega$.

\subsection{Discontinuous Galerkin discretization}
In this section we briefly discuss a Discontinuous Galerkin (dG) formulation associated with~\eqref{eq:weak_nonlin} which allows us to construct the adaptive stabilized formulation. We use as a starting point the dG discretization based on classical interior penalty schemes for elliptic problems~\cite{ern2006, cockburn2012}.
\subsubsection{Discrete setting}
Let $\{\mathscr{P}_h\}$ be a family of simplicial meshes of $\Omega$. For simplicity, we assume that $\Omega$ is exactly represented by any mesh in $\mathscr{P}_h$, that is, $\Omega$ is an interval, a polygon, or a polyhedron. We denote $T$ as the generic element in $\mathscr{P}_h$, with boundary $\partial T$, diameter $h_T$, and unit outward normal $\mathbf{n}_T$. We set $h = \max_{T \in \mathscr{P}_h} h_T$ and we assume, without loss of generality, that $h \leq 1$. We define the classical dG approximation space
\begin{align}
V_h :=\{v_h \in L^2(\Omega) \ | \  \forall T \in \mathscr{P}_h, v_h|_T \in \mathds{P}_k \},
\end{align}
where $\mathds{P}_k$ denotes the set of polynomials, defined over $T$, with polynomial degree smaller or equal than $k$. It is also convenient to set the extended space ${V_{h,\#} =H^2(\mathscr{P}_h)+V_h}$.

We collect all the faces or edges $F$ of $\Omega_h$ into the set ${\mathscr{S}_h= \bigcup_{T\in\mathscr{P}_h} F}$. We define the boundary skeleton $\mathscr{S}^\partial_h$ as ${\mathscr{S}^\partial_h = \mathscr{S}_h \cap \Gamma}$, and the internal skeleton $\mathscr{S}^0_h$ as ${\mathscr{S}^0_h=\mathscr{S}_h \backslash \Gamma}$.  Over $\mathscr{S}_h$, we define $\mathbf{n}_F$ as
a predefined normal over each $F$ being coincident with $\mathbf{n}$ when ${F \in \mathscr{S}_h^\partial}$, and $h_F$ as the diameter of the face $F$. On interior faces, for any function ${v \in V_{h, \#}}$, the jump  $[\![v]\!]_F$ and the standard (arithmetic) average $\{v\}_F$ are defined as
\begin{align*}
[\![v]\!]_F := v^- - v^+\quad \quad , \quad \quad \{v\}_F:= \frac{1}{2}(v^- + v^+),
\end{align*}
with $v^-$ and $v^+$ denoting the left and right face values respectively, with respect to the predefined normal $\mathbf{n}_F$. For ${F \in \mathscr{S}^\partial_h}$, we set ${[\![v]\!]_F=\{v\}_F=v|_F}$.

Finally, for a given norm $(\boldsymbol{\cdot},\boldsymbol{\cdot})_{V_h}$ of the discrete space $V_h$, we define the dual norm $\| \boldsymbol{\cdot} \|_{V^*_h}$ for $\phi \in V^*_h$ as:
\begin{align}\label{eq:dual-norm}
\| \phi \|_{V^*_h} := \sup_{v_h \in V_h \backslash \{0\}} \frac{\langle \phi, v_h \rangle_{V^*_h \times V_h}}{\| v_h \|_{V_h}},
\end{align}
where $\langle \boldsymbol{\cdot}, \boldsymbol{\cdot} \rangle_{V_h^* \times V_h}$ denotes the duality pairing in $V_h^* \times V_h$.

\subsubsection{Nonlinear discontinuous Galerkin formulation}
Considering the above discrete setting, we build the dG formulation for the continuous weak variational formulation  of~\eqref{eq:weak_nonlin} as
\begin{equation}\label{eq:dg_nonlin}
\left\{\begin{array}{l}
\text{Find } u_h \in V_h, \text{ such that:} \smallskip \\
\begin{array}{l}
n_h(u_h; v_h)= \ell_h(v_h), \quad \forall v_h \in V_h,
\end{array}
\end{array}\right.
\end{equation}
with
\begin{equation*}
\begin{array}{rl}
\displaystyle n_h(u_h;v_h) := & \displaystyle \sum_{T \in \mathscr{P}_h} (\nabla u_h \, , \, \nabla v_h)_{T} + \sum_{T \in \mathscr{P}_h} (\mathcal{F}(u_h)\, , \,v_h)_{T} \\
& + \displaystyle \sum_{F \in \mathscr{S}_h}\left[ \left( [\![ u_h ]\!] \, , \,  \{  \nabla v_h \} \cdot \mathbf{n}_F \right)_F - \left( \{  \nabla u_h \} \cdot \mathbf{n}_F  \, , \,  [\![ v_h ]\!]\right)_F +  \dfrac{\gamma}{h_F}\left( [\![ u_h ]\!], [\![ v_h ]\!] \right)_F \right], \\
\displaystyle \ell_h(v_h) :=& \displaystyle \sum_{T\in\mathscr{P}_h}(f, v_h)_T.
\end{array}
\end{equation*}
In the above, $\gamma>0$ is a user defined constant that we set as $\gamma = 3(k+1)(k+2)$, being $k$ the polynomial degree of the test space. Besides, we recall~\eqref{eq:der-def} and set the discrete derivative as:
\begin{equation}\label{eq:cnoidal_linearized}
\begin{array}{rcl}
\displaystyle n'_h(u_h; z_h, v_h)& :=&\displaystyle\frac{d}{d\epsilon}n_h(u_h+\epsilon z_h; v_h)\big|_{\epsilon=0} \vspace{0.2cm}\\
&=&\displaystyle \sum_{T \in \mathscr{P}_h} (\nabla z_h \, , \, \nabla v_h)_{T} + \sum_{T \in \mathscr{P}_h} (\mathcal{F}'(u_h) \, z_h\, , \,v_h)_{T} \\
& +&\displaystyle \displaystyle \sum_{F \in \mathscr{S}_h}\left[ \left( [\![ z_h ]\!] \, , \,  \{  \nabla v_h \} \cdot \mathbf{n}_F \right)_F - \left( \{  \nabla z_h \} \cdot \mathbf{n}_F  \, , \,  [\![ v_h ]\!]\right)_F +  \dfrac{\gamma}{h_F}\left( [\![ z_h ]\!], [\![ v_h ]\!] \right)_F \right],
\end{array}
\end{equation}
where
\begin{equation}
\mathcal{F}'(u_h)=\lambda \, m \, (1+u_h)^{m-1} - \mu \, \beta \exp(\beta u_h).
\end{equation}
We build our resolution scheme using~\eqref{eq:cnoidal_linearized}. This linearized form can be seen as a reaction-diffusion form in each increment $z_h$, for that reason, we provide the discrete space $V_h$ with a diffusion-type norm:
\begin{align}\label{eq:diff-norm}
\|w\|_{V_h}^2 :=& \; \theta \ \|w\|^2_{0, \Omega} + \| \nabla w\|^2_{0, \Omega} +  \sum_{F \in \mathscr{S}_h} \left(  \dfrac{\gamma}{h_F}[\![w]\!] , [\![w]\!] \right) _{0,F}.
\end{align}
with ${\theta=\lambda \, m \, A^{m-1}}$, where ${A>0}$ is a given constant associated to the maximum value of the normalized stress in the cnoidal solution.

\subsection{Adaptive stabilized finite element method based on residual minimization}
The stabilized finite element method based on residual minimization presented in~\cite{calo2019} delivers a mixed problem, with a saddle-point structure, in the case of linear problems. In this section, we develop the discrete formulation for the continuous problem~\eqref{eq:abs_vf}, which reads:
\begin{equation}\label{eq:nonlin}
{N_h(u_h)=\ell_h(\boldsymbol{\cdot})},
\end{equation}
where ${N_h : U_h \rightarrow V_h^*}$ represents the discrete nonlinear map with ${\langle N_h(z_h), v_h\rangle_{V_h^* \times V_h} := n_h(z_h; v_h)}$, being $n_h(\boldsymbol{\cdot}\, ;\boldsymbol{\cdot})$ the nonlinear form associated with problem~\eqref{eq:dg_nonlin}.

The methods solves for a \textit{continuous} approximation solution in a given discrete space $U_h$ (for instance, standard FEM functions), adequately considered as a subspace of the dG space $V_h$, i.e., ${u_h \in U_h = V_h \cap C^0(\Omega) \subset V_h}$. For example, in~\eqref{eq:dg_nonlin}, ${u_h \in U_h = V_h \cap H^1(\Omega)}$. The solution $u_h$ is then computed through minimizing the residual ${\ell_h(\cdot) - N_h (z_h)}$ associated to~\eqref{eq:nonlin}  in the norm of $V_h^*$:
\begin{align}\label{resmin_nonlin}
u_h = \operatorname*{argmin}_{z_h \in U_h} \frac{1}{2} \|\ell_h(\boldsymbol{\cdot}) - N_h (z_h)\|^2_{V_h^*} = \operatorname*{argmin}_{z_h \in U_h}\frac{1}{2} \|R_{V_h}^{-1} (\ell_h(\boldsymbol{\cdot})- N_h (z_h)) \|^2_{V_h},
\end{align}
with $\| \boldsymbol{\cdot} \|_{V_h}$ the norm defined in~\eqref{eq:diff-norm}, and $R^{-1}_{V_h}$ the inverse of the Riesz map (cf.,~\cite{oden2017}, Theorem 6.4.1). Similar to the thinking for~\cite{calo2019}, the nonlinear problem can be stated as a critical point of the minimizing functional, which translates into the following linear problem: Find $u_h \in U_h$ such that:
\begin{align}\label{eq:resmin_as_nonLinProb}
(R_{V_h}^{-1}(\ell_h - N_h(u_h)),R_{V_h}^{-1}DN_h(u_h; z_h))_{V_h} = 0, \quad \forall z_h \in U_h,
\end{align}
being $DN_h(u_h; z_h)$ the discrete form of the derivative (cf.,~\eqref{eq:der-def}). As noticed in~\cite{cohen2012}, problem~\eqref{eq:resmin_as_nonLinProb} can be equivalently written as the following saddle-point problem:

\begin{equation}\label{eq:nl-femwdg}
\left\{\begin{array}{l}
\text{Find } (\varepsilon_h, u_h) \in V_h \times U_h,  \text{ such that:} \smallskip \\
\begin{array}{rcll}
(\varepsilon_h, v_h)_{V_h} + n_h(u_h;v_h) &=& \ell_h(v_h), & \quad  \forall v_h \in V_h, \smallskip \\
n'_h(u_h; z_h, \varepsilon_h) &=& 0, &\quad  \forall z_h \in U_h.
\end{array}
\end{array}\right.
\end{equation}
The system~\eqref{eq:nl-femwdg} delivers simultaneously a stable and continuous approximation $u_h \in U_h$ of the dG formulation, and a residual representation $\varepsilon_h \in V_h$ that guides the adaptive mesh refinement. The first line of this system represents the nonlinear problem to solve, associated to the residual, whereas the second line can be seen as a constraint imposed on the tangent space built from the linearized form.

\subsubsection{Linearized problem}

We use Newton's method for solving the nonlinear problem. Given the discrete solution pair ${(\varepsilon^i_h, u^i_h)}$ of an iterative step $i$, we look for the increment ${(\delta \varepsilon_h, \delta u_h)}$ of the next iteration, and we set ${u}^{i+1}_h={u}^{i}_h+t^i {\delta u}_h$, and ${\varepsilon}^{i+1}_h={\varepsilon}^{i}_h+t^i {\delta \varepsilon}_h$, where $t^i$ represents a relaxation parameter to control the increment size. The method seeks for the solution pair $(\varepsilon_h^{i+1}, u_h^{i+1})$ that satisfies~\eqref{eq:nl-femwdg}. For the $i+1$-th iteration, the linearization of~\eqref{eq:nl-femwdg} reads:
\begin{equation}\label{eq:nl-newton}
\left\{\begin{array}{l}
\text{Given the pair } (\varepsilon^i_h, u^i_h), \text{ find } (\delta \varepsilon_h, \delta u_h) \in V_h \times U_h,  \text{ such that:} \smallskip \\
\begin{array}{rcll}
\hspace{-0.25cm}(\delta \varepsilon_h, v_h)_{V_h} + n'_h(u^i_h; \delta u_h, v_h) & \hspace{-0.25cm}= & \hspace{-0.25cm} \ell_h(v_h)- (\varepsilon^i_h, v_h)_{V_h} -n_h(u^i_h; v_h) , & \quad \forall v_h\in V_h,   \\
\hspace{-0.2cm}n'_h(u^i_h; z_h, \delta \varepsilon_h) & \hspace{-0.25cm}= & \hspace{-0.25cm}-n'_h(u^i_h; z_h, \varepsilon^i_h),                              & \quad \forall z_h \in U_h.
\end{array}
\end{array}\right.
\end{equation}

In matrix form, formulation~\eqref{eq:nl-newton} reads:
\begin{align}\label{eq:femwdg_mtx}
\begin{pmatrix}
\ G & B_u \ \\ \ B_u^T & 0 \
\end{pmatrix}
\begin{pmatrix}
\ \delta \varepsilon_h \ \\
\ \delta u_h \
\end{pmatrix}
=
\begin{pmatrix}
\ L  \ \\
\ 0 \
\end{pmatrix}
-
\begin{pmatrix}
\ G \varepsilon^i_h + N(u^i_h)\, \ \\
\ B^T_u \varepsilon^i_h \,\
\end{pmatrix}
\end{align}
where $G$ is the Grammian matrix built for the inner product that induces the norm in the discrete space $V_h$, $N(u^i)$ is the vector associated to the nonlinear form $n_h(u_h;v_h)$ and $B_u$ is the matrix associated with its linearization $n'_h(u_h^{i}; \delta u_h, v_h)$. The residual representative $\varepsilon_h$ is an implicit function of $u_h$. We define the pair ${\boldsymbol{ x_h}=( \varepsilon_h,  u_h)}$ that comprises both the solution and the residual representative, being valid also for the increments, which allows us to rewrite~\eqref{eq:femwdg_mtx} as:
\begin{equation*}
\boldsymbol{J}^i \, \boldsymbol{\delta x_h} = \boldsymbol{R}^i,
\end{equation*}
where
\begin{equation*}
\begin{array}{c}
\boldsymbol{J}^i =
\begin{pmatrix}
\ G & B_u \ \\ \ B_u^T & 0 \
\end{pmatrix}
\, \; \text{ and } \; \,
\boldsymbol{R}^i
=
\begin{pmatrix}
\ L  \ \\
\ 0 \
\end{pmatrix}
-
\begin{pmatrix}
\ G \varepsilon^i_h + N(u^i_h)\, \ \\
\ B^T_u \varepsilon^i_h \,\
\end{pmatrix}
\end{array}
\end{equation*}
The convergence of the method is controlled by the size of each iteration step through the relaxation parameter $t^i$. For that purpose, we use the damped Newton's method~\cite{bank1981} shown in Figure~\ref{fig:flow-chart}.

\begin{figure}[h!]
        \centering
        \includegraphics[width=0.8\textwidth]{simple-flow-chart.pdf}
        \caption{Damped Newton's algorithm flow chart.}
        \label{fig:flow-chart}
\end{figure}

\section{One-dimensional numerical experiments}\label{sec:numexp_1d}

In this section, we develop several one-dimensional numerical examples to illustrate the performance of the adaptive stabilized finite element method in the context of the cnoidal equation. Simulation results reported use FEniCS~\cite{ Alnaes2015} in \S \ref{ss:single-peak} \& \S \ref{ss:mult-peak}, and REDBACK~\cite{Tung2017} in \S \ref{ss:cont}. The main drawback of standard FEM implementations for this kind of problem lies in the difficulty to resolve the localized peak and particularly to find their location, which ends up delivering low-quality solutions. These limitations severely restrict its usage, forcing the initial guess to be close enough to the actual solution for the algorithm to converge, which is impractical. To overcome this these limitations in a nonlinear framework, we seek an algorithm with the ability to find automatically the location of the peaks. In practice, this means that we can start from an arbitrary initial trial solution with peaks located far from the final configuration. The numerical examples demonstrate that the adaptive stabilized finite element method (cf., \S\ref{sec:num}) can easily overcome these issues.

\subsection{Single peak solution}\label{ss:single-peak}

Thanks to the enhanced stability of the framework, we can solve~\eqref{eq:numsol_def_equation} and retrieve the expected peak solution for appropriate parameters. Figure~\ref{fig:comparison_mathematica} shows the comparison between the semi-analytical solution computed with Mathematica~\cite{wolfram1999mathematica} and the results obtained with our approach for ${\lambda=10,\, m=3, \,\mu=10^{-4}}$ and ${\beta=10}$, starting from an initial guess ${u_{IG}=2 \exp({-100(x-0.5)^2})}$ on a regular mesh of 100 nodes, getting to 273 nodes after four levels of adaptivity. We observe an excellent match, including the peak location, shape, and intensity, as shown in Figure~\ref{fig:1peak-zoom}.

\begin{figure}[h!]
    \centering
    \begin{subfigure}[c]{0.49\textwidth}
        \centering
        \includegraphics[width=\textwidth]{comparison_mathematica.pdf}
        \caption{\centering{Whole solution.}}
        \label{fig:1peak}
    \end{subfigure}
    \begin{subfigure}[c]{0.49\textwidth}
        \centering
        \includegraphics[width=\textwidth]{comparison_mathematica_zoom.pdf}
        \caption{\centering{Zoom around $x=0.5$.}}
        \label{fig:1peak-zoom}
    \end{subfigure}
    \caption{Comparison showing the good match between the Mathematica solution and our results for $\lambda=10$, starting from the initial guess ${u_{IG}=2 \exp(-100(x-0.5)^2)}$.}
    \label{fig:comparison_mathematica}
\end{figure}

\subsection{Multiple peak solution}\label{ss:mult-peak}
Given that we can recover known semi-analytic solutions with the method, we explore a harder problem considering more peaks. In this context, standard FEM approach is no longer appropriate, due to its lack of stability. For the following numerical examples, we define
\begin{equation}\label{eq:u_ig_1d}
{u_{IG}(x) := A_0 \, \left[\frac{\exp({-1250(x-x_0)^2})}{\sin(x_0\pi)}+ \frac{\exp({-1250(x-(1-x_0))^2})}{\sin((1-x_0)\pi)}\right]\sin(x\pi) }
\end{equation}
as the initial guess function, being $x_0$ the arbitrary location of the first peak ${(0 \leq x_0\leq 0.5)}$. This initial guess choice implies that the second peak location is at $1-x_0$. Besides, $A_0>0$ is an arbitrary number that coincides with the values of $u$ at the peak locations in the initial guess. Table~\ref{tab:param} shows the set parameters for the examples in this subsection.

\begin{table}[ht!]
	\small
	\centering
	\caption{Parameters for the one-dimensional numerical examples with two-peak solution.}
	\begin{tabular}{|c|c|c|c|c|c|c|}
		\hline
	 	Example&$\lambda$ &$m$ & $\mu$& $\beta$ & $A_0$ & $x_0$ \\
	 	\hline
		3.2.1&40&3&$10^{-4}$&10&1.80&0.200\\
		\hline
		3.2.2&40&3&$10^{-4}$&10&1.80&0.175\\
		\hline
		3.2.3&40&3&$10^{-4}$&10&1.80&0.425\\
		\hline
		3.2.4&40&$\pi$&$10^{-4}$&10&2.70&0.350\\
		\hline
	\end{tabular}
	\label{tab:param}
\end{table}

\subsubsection{Comparison against standard FEM}
As a first example, Figure~\ref{fig:sol-2peaks} shows the results comparing the discrete solution obtained with the standard FEM formulation and the new adaptive stabilized method, for the same initial guess. For this two-peak example, we set the arbitrary location as ${x_0=0.2}$ (see Table~\ref{tab:param}). We use cubic trial functions $(\mathds{P}_3)$ for both methods, but we take advantage of the possibility of enriching the test space in the adaptive stabilized method, using test functions one degree higher $(\mathds{P}_4)$. Finally, we use a fixed mesh for the standard finite element solution of mesh size of ${h = 10^{-6}}$, to develop a fair comparison with the final refined mesh obtained through the adaptive method, which starts from a mesh size of 100 elements (${h=0.01}$) and gets no finer than ${h=10^{-6}}$ locally. We can observe that this new technique properly captures the final location of the peaks at $x\approx\{0.27, 0.73\}$, whereas the standard method delivers spurious oscillations.
\begin{figure}[h!]
	\centering
	\includegraphics[width=0.6\textwidth]{solution_2peaks.pdf}
	\caption{Comparison between standard finite element solution (in red) and the adaptive stabilized method (in green) for an initial guess ${u_{IG} := 1.8 \, [\exp({-1250(x-0.2)^2})/\sin(0.2\pi)+ \exp({-1250(x-0.8)^2})/\sin(0.8\pi)]\sin(x\pi)}$ (in blue). The stabilized method converges to a stable solution, whereas the standard FEM approach leads to spurious oscillations.}
	\label{fig:sol-2peaks}
\end{figure}

\subsubsection{Initial guess close to the boundaries}
We now investigate examples with different conditions to show that the new method converges robustly with respect to the initial condition. The distance  between the initial and final peak locations using the adaptive method is significantly larger than that of the one for which standard FEM on a fine mesh can converge. In this example, we locate the peaks close to the boundaries (see Table~\ref{tab:param}). Figure~\ref{fig:2peaks-175} shows the iterative solutions obtained at each refinement step. The adaptive method converges when standard FEM fails, even with an order of magnitude finer mesh (${h=10^{-7}}$). The adaptive approach successively corrects the peak locations at each refinement level, starting from a mesh size of $h=0.01$. After 47 refinement levels (approximately, $24,000$ iterations), we obtain a solution with a final residual norm $\|\boldsymbol{R}^{i+1}\|_0<10^{-9}$. Incidentally, this example also shows that the solution can be asymmetric.
\begin{figure}[h!]
	\centering
	\includegraphics[width=0.6\textwidth]{x=0175.pdf}
	\caption{Evolution of the profiles at intermediate mesh refinement steps for an initial guess using ${x_0=0.175}$ (in blue) showing convergence to an appropriate (asymmetrical) solution (in red).}
	\label{fig:2peaks-175}
\end{figure}

\subsubsection{Initial guess close to the center}
In this example, we locate the initial peaks close to the center (see Table~\ref{tab:param}). In that instance, we end up converging to the solution of Figure~\ref{fig:2peaks-425} using the adaptive stabilized method after 37 refinement levels (approximately, $7,400$ iterations) for the same tolerance than in the previous case. The final peak locations are $x\approx\{0.27, 0.73\}$. FEM is not able to converge using this initial guess either. From our experience, which we do not report for the sake of brevity, we find that FEM simulations require initial guesses sufficiently close to the final solution for the method to converge. In practice, FEM requires the distance between the initial and the final solution peak to be at least an order of magnitude smaller than what the adaptive stabilized method admits.
\begin{figure}[h!]
	\centering
	\includegraphics[width=0.6\textwidth]{x=0425.pdf}
	\caption{Evolution of the profiles at intermediate mesh refinement steps for an initial guess using ${x_0=0.425}$ (in blue), showing convergence to an appropriate (symmetrical) solution (in red).}
	\label{fig:2peaks-425}
\end{figure}

\subsubsection{Non-integer exponent}
Finally, we simulate a scenario with a non-integer exponent $m=\pi$ to show the robustness of the adaptive stabilized method for an irrational exponent. As mentioned in \S\ref{sec:num}, the analytical approach of the cnoidal equation does not provide solutions for this class of exponents. We consider the initial guess~\eqref{eq:u_ig_1d} with $x_0=0.35$ and we set higher peaks values than in past examples (see Table~\ref{tab:param}). Figure~\ref{fig:sol-pi} displays the evolution of the solution profile, showing that the method can robustly simulate irrational exponents larger than 3, a limitation of the analytical resolution approach~\cite{Veveakis2015}. Also, the adaptive stabilized scheme correct the height of the peak values. Numerical simulations not presented here also showed a good performance for even higher exponent values up to 7.
\begin{figure}[h!]
    \centering
    \includegraphics[width=0.6\textwidth]{m=pi.pdf}
    \caption{Evolution of the two-peak solution profile using the adaptive stabilized method for $m=\pi$ and an initial guess with peak values of ${u(0.35)=u(0.65)=2.7}$ (in blue), showing convergence to an appropriate (symmetrical) solution (in red).}
    \label{fig:sol-pi}
\end{figure}

\subsection{Advantages of the usage of the adaptive FEM method in the cnoidal problem}
In summary, the adaptive stabilized method can robustly and efficiently predict the solution patterns for many relevant configurations for the cnoidal wave problem. This heretofore unsurmountable problem is now solvable and the performance of the method can be explained for two reasons. Firstly, this technique allows the enrichment of the test space, which implies having stronger approximation norms and better behavior in the solution at each level. This enrichment is not possible for FEM, due to the structure of Galerkin's approach. Secondly, the adaptive mesh refinement scheme ensures the elimination of oscillations through the reduction of the local error. This robust adaptivity represents an important feature of the method, especially because of the localized nature of the solution. Using this adaptive stabilized method, we can build an indicator from the residual representative that allows us to correct the location of the peaks at each refinement level.

\subsection{Numerical continuation}\label{ss:cont}

We can now investigate in detail the influence of the parameter $\lambda$, to better understand the most notable aspects of the solution. In particular, we are interested in identifying the range of $\lambda$ values where a peak solution exists. More precisely, we seek to discern the presence of a threshold $\lambda_c$, above which multiple peaks can be obtained, as well as the relationship between the number, $p$, of peaks and $\lambda$. For example,~\cite{Regenauer-Lieb2013a} reported a ${p\propto\sqrt{\lambda}}$. To that effect, we use continuation technique to study the stability regimes of the system, starting from identified solutions.

We note that it is possible to solve the cnoidal equation~\eqref{eq:cnoidal_1} in the stationary limit via standard FEM procedures starting from a sufficiently accurate initial solution, which allows us to retrieve the continuation map using classical techniques. Hence, we perform stability analysis using a pseudo-arclength continuation algorithm from~\cite{Keller1977}, which allows us to identify different behavior characteristics of the solution of the system (\ref{eq:cnoidal_chem}) where we also consider $\lambda$ as a variable. The algorithm starts with two solutions $(u_{0}, u_1)$ of (\ref{eq:cnoidal_chem}) for the respective values $(\lambda_0, \lambda_1)$, where each solution is computed on a discretized mesh of $n$ points. The next step is then computed starting from the initial guess ${[2u_1-u_{0},2\lambda_1-\lambda_0]^T}$, i.e. following the tangent vector of $[u,\lambda]^T$ in $\mathbb{R}^{n+1}$ and searching for a solution $[u,\lambda]^T$ at a distance $\Delta s$ from the previous one. We compute following steps recursively. The extended system of equations solved can then be written as
\begin{subequations}
    \label{eq:continuation_system}
    \begin{align}
    & \frac{\partial^2 u_n}{\partial z^2}-\lambda(1+u_n)^{m}+\mu \mathrm{e}^{\beta u_n}=0, \\
    &\dot{u}_{0}\left(u-u_{0}\right)+\dot{\lambda}_{0}\left(\lambda-\lambda_{0}\right)=\Delta s,
    \end{align}
\end{subequations}
where $[\dot{u}_{0},\dot{\lambda}_{0}]^T$ is the tangent vector at $[u_{0},\lambda_{0}]^T$.

\subsubsection{One peak solutions}

Solving (\ref{eq:cnoidal_chem}) for any positive value of $\lambda$ with a zero initial guess leads to the traditional Terzaghi's consolidation isochrone profile $u(x)\leq0,  \: \forall x \in [0,1]$~\cite{Terzaghi1943} shown in Figure~\ref{fig:ccurve-point-b}. Running the numerical continuation algorithm from two such initial solutions for increasing values of $\lambda$ yields the expected branch of traditional Terzaghi's solutions, with the minimum value $u_c$, obtained at the center ($x=0.5$) for symmetry reasons, decreasing asymptotically towards $-1$ for ${\lambda \to \infty}$. Running the continuation algorithm for decreasing values of $\lambda$ provides the evolution of the solution as $\lambda$ becomes negative. Figure~\ref{fig:c-curve-plot} shows the results plotted in the ${\lambda-u_c}$ space. We observe a C-shape solution, marked by a threshold value of $\lambda^{(1)}_{\min}\approx -1.42$ and the existence of an upper branch. The profiles of two solutions obtained for the same value of $\lambda=1.5$  but on the different branches, upper and lower, are also shown in Figs.~\ref{fig:ccurve-point-b}~and~\ref{fig:ccurve-point-c} respectively.
\begin{figure}[t!]
    \centering
    \begin{subfigure}[c]{0.8\textwidth}
        \centering
        \includegraphics[width=\textwidth]{C_curve_1_peak_A.pdf}
        \caption{}
        \label{fig:c-curve-plot}
    \end{subfigure}
    \begin{subfigure}[c]{0.4\textwidth}
        \vspace{0.5cm}
        \centering
        \includegraphics[width=\textwidth]{profile_lower.pdf}
        \caption{}
        \label{fig:ccurve-point-b}
    \end{subfigure}
    \begin{subfigure}[c]{0.4\textwidth}
        \vspace{0.5cm}
        \centering
        \includegraphics[width=\textwidth]{profile_upper.pdf}
        \caption{}
        \label{fig:ccurve-point-c}
    \end{subfigure}
    \caption{(a) C-curve result of the stability analysis of (\ref{eq:cnoidal_chem}), plotting the value $u_c$ of the solution at the center of the domain ($x=0.5$) with respect to parameter $\lambda$, capturing the traditional Terzaghi's consolidation (blue lower branch) and single peak case (red upper branch). Points B and C mark two points on the lower and upper branches for the same value of $\lambda=1.5$ and their respective solution profiles are shown in subfigures (b) and (c); (b) Terzaghi's solution profile for $\lambda=1.5$; (c) one-peak solution profile for $\lambda=1.5$}
    \label{fig:C-curve-1peak}
\end{figure}

\subsubsection{Multiple peak solutions}

The numerical continuation can also be applied in other areas of the solution space to trace different branches, whose presence is indicated by previous studies~\cite{Veveakis2014,Regenauer-Lieb2013a} who noted the existence of several peaks in the solution for higher values of $\lambda$. For the sake of demonstrating the existence and behavior of those solutions, we arbitrarily limit ourselves to solutions with up to seven peaks for ${\lambda<500}$.

For all given numbers of peaks ${(2 \leq  p \leq 7)}$, we apply the same procedure and identify manually a couple of solutions for values of $\lambda$ close to $500$ before running the continuation algorithm. Our procedure may not capture of all possible branches in the parameter space. Therefore, potentially even richer solutions exist. While the peak is always located at the center for single-peak solutions, as shown in Figure~\ref{fig:C-curve-1peak}, locations of the peaks do vary with $\lambda$ when there exists more than one peak.

\begin{figure}[h!]
    \centering
    \includegraphics[width=0.9\textwidth]{space-map.pdf}
    \caption{Map of the seven continuation curves obtained for $\lambda \leq 500$. All curves, for $p$ peaks, display a C-shape, marked by a minimum value $\lambda^{(p)}_{\min}$, as shown in Figure~\ref{fig:C-curve-1peak} for $p=1$.}
    \label{fig:space_map}
\end{figure}

Figure~\ref{fig:space_map} shows the corresponding results obtained for $\lambda \leq 500$, plotting the solution amplitude for each point of the curves. Similarly to the single-peak branch ($p=1$), all branches for $p\geq 2$ display a C-shape, indicating the existence of a lower threshold $\lambda^{(p)}_{\min}$ limiting the validity region of each branch. Each C-shape curve also indicates that there are at least two possible profiles for any given number $p$ of peaks for any admissible value of $\lambda \geq \lambda^{(p)}_{\min}$, which differ by the shape of the solution in general, with the same maximum value but different spacings between peaks.

\begin{figure}[ht!]
    \centering
    \begin{subfigure}[c]{0.75\textwidth}
        \centering
        \includegraphics[width=\textwidth]{peak_spacing.pdf}
        \caption{Spacing between peaks for branches with $p$ peaks (${2\leq p \leq7}$), with two values of $\lambda$ highlighted: ${\lambda^{(2)}_{\min}}$ (marking the lowest value of $\lambda$ on the C-curve) and ${\lambda^{(2)}_\text{eqd}}$ (value of $\lambda$ at point C, where ${spacing=1/p}$).}
        \label{fig:spacing-a}
    \end{subfigure}
    \begin{subfigure}[c]{0.26\textwidth}
        \centering
        \includegraphics[width=\textwidth]{spacing_point_B.pdf}
        \caption{Point B: far-peaks sol.}
        \label{fig:spacing-b}
    \end{subfigure}
    \begin{subfigure}[c]{0.26\textwidth}
        \centering
        \includegraphics[width=\textwidth]{spacing_point_C.pdf}
        \caption{Point C: equidistant-peaks sol.}
        \label{fig:spacing-c}
    \end{subfigure}
    \begin{subfigure}[c]{0.26\textwidth}
        \centering
        \includegraphics[width=\textwidth]{spacing_point_D.pdf}
        \caption{Point D: close-peaks sol.}
        \label{fig:spacing-d}
    \end{subfigure}
    \begin{subfigure}[c]{0.75\textwidth}
        \centering
        \includegraphics[width=\textwidth]{inv_peak_spacing.pdf}
        \caption{Inverse spacing plot, highlighting for each C-curve with $p$ peaks the turning point  (marked by point C for ${p=2}$) where ${1/spacing=p}$,  as well as the asymptotic behavior ${\lim_{\lambda\to\infty} spacing=p+1}$ on the upper branch (point B) and  ${\lim_{\lambda\to\infty} spacing=p+1}$ on the lower branch (point D). }
        \label{fig:spacing-e}
    \end{subfigure}
        \caption{Spacing maps as a function of $\lambda$ for every branch of the $p$-peak solution (${2 \leq p \leq  7}$).}
        \label{fig:peak_spacing}
    \end{figure}

We compute the distance between peaks by post-processing, noting that the regularity of all solutions provides constant spacing for the various peaks on any given profile. Figure~\ref{fig:peak_spacing} shows the spacing evolution (and its inverse) as a function of $\lambda$ for every branch of the $p$-peak solution for $2 \leq p \leq  7$, with all branches following the same pattern. For example, the two-peak solution on Figures~\ref{fig:spacing-b}-\ref{fig:spacing-d}. Figure~\ref{fig:spacing-a} shows the real spacing, Figure~\ref{fig:spacing-e} the inverse of that spacing, appears in Figures~\ref{fig:spacing-b}-\ref{fig:spacing-d} three specific profiles at points B, C and D highlighted on Figures~\ref{fig:spacing-a}~\&~\ref{fig:spacing-e}.

The equidistant spacing of $1/p$  (for ${x \in [0,1]}$) for all $p$ peaks is obtained for ${\lambda=\lambda^{(p)}_\text{eqd}}$, slightly larger than ${\lambda^{(p)}_{\min}}$, with the corresponding point on the C-curve for $p=2$ highlighted as point C on Figures~\ref{fig:spacing-a}~\&~\ref{fig:spacing-e}. Figure~\ref{fig:spacing-c} shows the corresponding solution profile. This point C naturally separates the higher and lower branches, leading to point B (resp. D) with increasing (resp. decreasing) spacing. Figures~\ref{fig:spacing-e} shows the asymptotic behavior ${\lim_{\lambda\to\infty} \text{spacing}=p+1}$ on the upper branch (point B) and  ${\lim_{\lambda\to\infty} \text{spacing}=p+1}$ on the lower branch (point D). The three solution profiles at points B, C and D are shown in Figs.~\ref{fig:spacing-b}-\ref{fig:spacing-d} and highlight this spacing evolution across the two-peak branches. All other branches (for ${p>2}$) in Figures~\ref{fig:peak_spacing} display the same behavior.

\subsubsection{Evolution of number of peaks with respect to $\lambda$}
The identification of the natural turning point for ${\lambda=\lambda^{(p)}_\text{eqd}}$  on the lower branch of each C-curve, which is slightly different from the point marking the minimum possible of ${\lambda}$ ${(\lambda^{(p)}_\text{min})}$, indicates that both characteristic points are useful to refer to the ``start" of a C-curve. Table~\ref{tab:lambdas} lists the numerical values, with a precision of two decimals, of ${\lambda^{(p)}_\text{min}}$ and ${\lambda^{(p)}_\text{eqd}}$ for ${1\leq p \leq 7}$ as shown on Figure~\ref{fig:space_map}.
Figure~\ref{fig:sqrt_lambda} plots those evolutions of ${\lambda=\lambda^{(p)}_\text{min}}$ and ${\lambda=\lambda^{(p)}_\text{eqd}}$ with respect to $\lambda$. Fitting of those curves shows correspondence for both cases with the square-root relationship previously identified in~\cite{Regenauer-Lieb2013a} between $\lambda$ and the number of peaks $p$.

\begin{table}[ht!]
    \small
    \centering
    \caption{Values for each branch with $p$ peaks in Figure~\ref{fig:peak_spacing}, of ${\lambda^{(p)}_{\min}}$, the minimum value of $\lambda$ reached, and ${\lambda^{(p)}_\text{eqd}}$, the value of $\lambda$ for which the peaks are equidistant with ${spacing=1/p}$ (naturally, no such value exists for ${p=1}$).}
    \begin{tabular}{|c|c|c|c|c|c|c|c|}
        \hline
        $p$ &1 & 2 & 3 & 4 & 5 & 6 & 7 \\
        \hline
        $\lambda^{(p)}_\text{min}$ &-1.42&26.99& 66.21&121.93&194.38&283.70&390.00\\
        $\lambda^{(p)}_\text{eqd}$ &--&30.30&69.48&125.32&197.84&287.23&393.59\\
        \hline
    \end{tabular}
    \label{tab:lambdas}
\end{table}

\begin{figure}[ht!]
    \centering
    \includegraphics[width=0.475\textwidth]{sqrt_lambda.pdf}
    \includegraphics[width=0.475\textwidth]{sqrt_lambda_eqd.pdf}
    \caption{Evolution of $p$, the number of peaks, with respect to $\lambda^{(p)}_{\min}$ and  $\lambda^{(p)}_\text{eqd}$ as defined in Figure~\ref{fig:peak_spacing}.}
    \label{fig:sqrt_lambda}
\end{figure}

\section{Discussion and conclusions}
\label{sec:discussion}

Regarding the chemical effects in the physical formulation, following the work from~\cite{Alevizos2017} provides a bounded solution of the original problem of cnoidal waves in solids~\cite{Veveakis2015,Regenauer-Lieb2013a}. We solve the resulting nonlinear equation using an adaptive stabilized finite element framework~\cite{calo2019} for a wide spectrum of scenarios. We adopt a numerical continuation analysis to map the influence of the diffusivity ratio $\lambda$ in (\ref{eq:cnoidal_chem}) on the solution behavior, providing new insights on this recently described volumetric instability. Figure~\ref{fig:C-curve-1peak} shows that we can retrieve both the classical Terzaghi's isochrone (Figure~\ref{fig:ccurve-point-b}) or a peak-type solution (Figure~\ref{fig:ccurve-point-c}) for positive values of $\lambda$ and not only inside a specific range of integer numbers, as implied in~\cite{Veveakis2015,Regenauer-Lieb2013a}. This outcome represents a starting point for further experimental investigations to check whether this response can be observed in laboratory tests.

We also recover the previously identified square root relationship of the number of peaks in the solution with $\lambda$~\cite{Veveakis2015,Regenauer-Lieb2013a} and analyze in detail the solution space, (cf., see Figure~\ref{fig:space_map}). The C-shape nature of all $p$-peak solutions identifies minimum thresholds ($\lambda^{(p)}_{\min}$) above which those solutions exist and Figure~\ref{fig:sqrt_lambda} shows that $\lambda^{(p)}_{\min}\propto\sqrt{\lambda}$. Figure~\ref{fig:space_map}, however, also shows that all $p$-peak solutions exist in a stable manner past the $\lambda$ value of $\lambda^{(p+1)}_{\min}$, where the $(p+1)$-peak solution appears. This result has important consequences for the physical understanding of the problem as the number of bands measured on natural observations cannot be automatically correlated with $\lambda$ as~\cite{Veveakis2015} states. Furthermore, Figure~\ref{fig:peak_spacing} shows that the spacing between bands follows itself a more complex pattern than initially thought, since a $p$-peak solution exists for all values of $\lambda\geq\lambda^{(p)}_{\min}$ on two branches (see Figure~\ref{fig:space_map}), with the spacing between peaks increasing with $\lambda$ (asymptotically towards $1/(p-1)$) on the lower branch and decreasing with $\lambda$ (asymptotically towards $1/(p+1)$) on the upper branch. Figure~\ref{fig:peak_spacing} also shows that the natural turning point for each C-curve with $p$ peaks corresponds to the solution profile with all peaks equidistant at a spacing of $1/p$, obtained for $\lambda=\lambda^{(p)}_\text{eqd}$. The number of peaks $p$ follows as well a square root relationship concerning those ($\lambda^{(p)}_\text{eqd}$) values.

All simulations assume Dirichlet boundary conditions of (\ref{eq:cnoidal_chem}). Thus, the most appropriate boundary conditions that can match the natural observations are still open~\cite{Veveakis2015,Veveakis2014}. We deploy more advanced bifurcation analysis tools to complete the bifurcation map and identify all branches in the parameter space, including unstable ones, along with their characteristic points. Nonetheless, the surprisingly rich information obtained through introducing an appropriate numerical scheme to a simple generalization of the theory of consolidation, suggests that future works using this scheme in more elaborate elasto-viscoplastic formulations could enhance the mechanical solution with additional modes of localization stemming from the volumetric part of the plastic increment.

\section*{Acknowledgments}
This research was partially supported by the Australian Government through the Australian Research Council's Discovery Projects funding scheme (projects  DP170104550, DP170104557). This publication was also made possible in part by the CSIRO Professorial Chair in Computational Geoscience at Curtin University and the Deep Earth Imaging Enterprise Future Science Platforms of CSIRO. Additional support was provided by the European Union's Horizon 2020 Research and Innovation Program of the Marie Skłodowska-Curie grant agreement No. 777778, and the Mega-grant of the Russian Federation Government (N 14.Y26.31.0013). At Curtin University, The Institute for Geoscience Research (TIGeR) and by the Curtin Institute for Computation, kindly provide continuing support. MV acknowledges support by the DE-NE0008746-DoE project.