# Mathematical Model for Transmission Dynamics of Tuberculosis in Burundi

Steve Sibomana  Kelly Joelle Gatore Sinigirira<EMAIL_ADDRESS>Paterne
Gahungu  David Niyukuri  Department of Mathematics, Faculty of Science,
University of Burundi, Bujumbura, Burundi  Doctoral School, University of
Burundi, Bujumbura, Burundi  Institute of Applied Statistics, University of
Burundi, Bubanza, Burundi  African Institute for Mathematical Sciences (AIMS),
Rwanda  The South African Department of Science and Technology–National
Research Foundation (DST-NRF) Centre of Excellence in Epidemiological
Modelling and Analysis (SACEMA), Stellenbosch University, Cape Town, South
Africa

###### Abstract

Tuberculosis (TB) is among the main public health challenges in Burundi. The
literature lacks mathematical models for key parameter estimates of TB
transmission dynamics in Burundi. In this paper, the supectible-exposed-
infected-recovered (SEIR) model is used to investigate the transmission
dynamics of tuberculosis in Burundi. Using the next generation method, we
calculated the basic reproduction number, $R_{0}$. The model is demonstrated
to have a disease-free equilibrium (DEF) that is locally and globally
asymptotically stable. When the corresponding reproduction threshold quantity
approaches unity, the model enters an endemic equilibrium (EE). That means,
the disease can be controlled through different interventions in Burundi. A
sensitivity analysis of the model parameters was also investigated. It shows
that the progression rate from latent to becoming infectious had the highest
positive sensitivity, which means that $R_{0}$ increases and decreases
proportionally with an increase and a decrease of that progression rate.

###### keywords:

Tuberculosis, Reproduction number, disease free-equilibrium, endemic
equilibrium, Lyapunov function, Burundi

††journal: Epidemics

##  1 Introduction

Tuberculosis (TB) is an airborne disease caused by Mycobacterium tuberculosis
(MTB)[1], and it is reported to be the second leading cause of morbidity and
mortality in the world from a single infectious agent after the human
immunodeficiency virus (HIV)[2]. The MTB spreads through inhaling droplets
from the cough or sneeze of a person suffering from active tuberculosis. The
bacteria enters the body causing an MTB infection affecting mainly the lungs,
but it can also affect any other part of the body including the urinary tract,
brain, lymph nodes, bones, joints and the ear. Person with lowered immunity
such as those with HIV, diabetes, immune disorders, and stage renal disease,
those on drugs that suppress immunity, young children, pregnant women among
others are at a higher risk of contracting the disease[3, 4].

In Burundi, TB remains a public health problem, it is rampant in an endemo-
epidemic mode and all layers of the population are concerned. In general,
infection with tuberculosis is very likely to be asymptomatic for healthy
people. The lifetime risk of developing clinically active TB after being
infected is about $10\%$ [5]. People who have latent TB infection are not
clinically ill or capable of transmitting TB [5, 6]. The immunity of older
people who have previously been infected may decrease, and they may then be at
risk of developing active TB of either exogenous reinfection (that means
acquisition of a new infection from another infectious individual) or an
endogenous reactivation of latent bacilli (that means the reactivation of a
pre-existing dormant bacillus infection)[7]. Latent and active tuberculosis
can be treated with antibiotics. However, its treatment has a side effects
(sometimes quite serious) and takes a long time.

Carriers of the tuberculosis bacillus who have not developed tuberculosis
disease can be treated with only one Isoniazid, also known as isonicotinic
acid hydrazide (INH). For active tuberculosis it is often used together with
rifampicin, pyrazinamide, and either streptomycin or ethambutol[8].
Unfortunately, it should be taken religiously for 6-9 months. Treatment of
people with active tuberculosis requires simultaneous use of three drugs for a
period of at least 12 months. Lack of compliance with these drug treatments,
do not only can cause a relapse, but the development of antibiotic resistant
tuberculosis, one of the most serious public health problems facing today’s
society.

Studying the spread of TB using statistics and mathematics models did not
receive enough attention in Burundi. As a result, we only observed a very
limited use of mathematical models in studying the dynamics of TB transmission
in Burundi’s human populations. In the broad scientific literature,
communicable diseases such as measles, influenza, rubella, among others, have
been studied by many mathematical models [9, 10, 11]. These diseases have a
number of traits, such as the fact that they frequently create epidemics and
that transmission rates are greatly influenced by age-dependent contact rates.
The etiological agents of these communicable diseases are viruses of different
families but all capable of generating similar symptoms. Waaler and Anderson
were the first who modeled mathematically tuberculosis transmission
dynamics[12].

In this work, a mathematical model for the transmission dynamics of TB in
Burundi has been developed. We determine the existence and positivity of the
system, and we provide the disease’s equilibrium points and the reproduction
number. The stability of the disease equilibrium is then determined. A
sensitivity analysis is performed on the model parameters. Based on the TB
data from Burundi, simulations and interpretations of the results will be
carried out.

##  2 Methods

###  2.1 Mathematical model

A mathematical model is be established using Susceptible-Exposed-Infected-
Recovered (SEIR) compartmental approach as shown on Figure 1, where $S(t)$,
susceptible humans, $E(t)$, exposed humans to tuberculosis, $I(t)$, infected
humans with active tuberculosis, $R(t)$, recovered humans at time $t\geq 0$.
We consider that susceptible population keep increasing by the human birth at
$h$ rate. Thee loose of immunity is considered to be at rate $\rho$, and there
is no permanent immunity to tuberculosis. Humans can contract MTB tuberculosis
through contact with individuals who are infected with the disease. Therefore,
they enter the exposed (latent) with $\beta$ rate. A proportion of the exposed
class develop active tuberculosis, thus, moving into the infectious class with
$\varepsilon$ rate of progression. If treatment is administered promptly,
those who recover from the disease will move to the recovered class at
$\gamma$ rate. Each human class decreases also by natural mortality at $\mu$
rate, except the infectious class which adds the mortality caused by the MTB
tuberculosis at $\tau$ rate. $N(t)$ is used to denote the total population at
time $t$ such that $N(t)=S(t)+E(t)+I(t)+R(t)$. All parameters of the model are
assumed to be non-negative and the total human is assumed to be constant.

Figure 1: Diagram of the Susceptible-Exposed-Infected-Recovered (SEIR) Model

From Figure 1, we can write the following system of ordinary differential
equations:

$\left\\{\begin{array}[]{ccc}\frac{dS(t)}{dt}&=&hN(t)-\mu S(t)-\beta
S(t)I(t)+\rho R(t),\\\ &&\\\ \frac{dE(t)}{dt}&=&\beta
S(t)I(t)-(\mu+\varepsilon)E(t),\\\ &&\\\ \frac{dI(t)}{dt}&=&\varepsilon
E(t)-(\mu+\tau+\gamma)I(t),\\\ &&\\\ \frac{dR(t)}{dt}&=&\gamma
I(t)-(\mu+\rho)R(t).\end{array}\right.$ (1)

The description of parameters and variables is the following

  1. Variables

    1. $S(t)$: susceptible humans,

    2. $E(t)$: exposed humans to tuberculosis,

    3. $I(t)$: infected humans with active tuberculosis,

    4. $R(t)$: recovered humans,

    5. $N(t)$: total population.

  2. Parameters

    1. $h$: human birth rate

    2. $\beta$: rate at which the susceptibles become exposed to MTB

    3. $\rho$: progression rate from recovered to susceptible class

    4. $\varepsilon$: progression rate from latent to infectious class

    5. $\gamma$: progression rate from recovered to susceptible class

    6. $\mu$: natural human death rate

    7. $\tau$: human death rate caused by the MTB tuberculosis

###  2.2 Existence and Positivity of Solutions

The system 1 is epidemiologically and mathematically well-posed in the
feasible region $\Gamma$ given by

$\Gamma=\\{(S,L,I,R)\in\mathbb{R}_{+}^{4}:~{}0<\mu\leq h\\}.$ (2)

######  Proposition 2.1.

The system (1) always admits positive solutions for all positive initial
conditions and the biological region $\Gamma\in\mathbb{R}_{+}^{4}$ is
positively invariant and globally attractive for the same system .

Proof:

From the system 1, we have

$\left\\{\begin{array}[]{ccl}\frac{dS(t)}{dt}&\geq&-(\mu+\beta I(t))S(t),\\\
&&\\\ \frac{dE(t)}{dt}&\geq&-(\mu+\varepsilon)E(t),\\\ &&\\\
\frac{dI(t)}{dt}&\geq&-(\mu+\tau+\gamma)I(t),\\\ &&\\\
\frac{dR(t)}{dt}&\geq&-(\mu+\rho)R(t).\end{array}\right.$ (3)

Integrating each inequality of the system 3 yields

$\left\\{\begin{array}[]{ccl}S(t)&\geq&S(0)\exp\\{-(\mu
t+\beta\int_{0}^{n}I(n)dn)\\},\\\ &&\\\
E(t)&\geq&E(0)\exp\\{-(\mu+\varepsilon)t\\},\\\ &&\\\
I(t)&\geq&I(0)\exp\\{-(\mu+\tau+\gamma)t\\},\\\ &&\\\
R(t)&\geq&R(0)\exp\\{-(\mu+\rho)t\\}.\end{array}\right.$ (4)

where $S(0),E(0),I(0),R(0)$ are all positive initial conditions. From the
system (4), all state variables [that means ($S(t),E(t),I(t),R(t)$) ] are all
positive for all $t\geq 0$. Thus, the solutions of model 1 remain positive in
$\Gamma$ for all time $t\geq 0$.

Add member to member the equations of the system (1), we notice that the total
of the human population satisfies the following relation:

$\dot{N}(t)=hN(t)-\mu N(t)-\tau I(t).$ (5)

We have $\dot{N}(t)\leq hN(t)-\mu N(t)$ according the equation 5. Then,
$\dot{N}(t)\leq 0$ if and only if $h\leq\mu$. So, $N(t)\leq
N(0)\exp[(h-\mu)t]$. We conclude that the region $\Gamma$ is positively
invariant. Moreover, if $h<\mu$, we see that the solution of the system 1
enter in $\Gamma$ in finite time, that means $N(t)$ asymptotically approaches
$h$.

Therefore, all solutions in $\mathbb{R}_{+}^{4}$ eventually enter $\Gamma$
that is the biological region $\Gamma$ is globally attractive for the same
system. So the problem is then mathematically and epidemiologically well
posed. Hence, every solution of the model 1 with initial conditions in
$\Gamma$ remains in $\Gamma$ for all $t\geq 0$. $\Box$

###  2.3 Disease Equilibria points

####  2.3.1 Disease Free Equilibrium (DFE)

A disease free equilibrium point is a solution of the system (1) in holding
that there is no disease in population. In this case $E=I=R=0$. Therefore in
our case the DFE, $E_{0}$, is expressed as

$E_{0}=\left(\frac{Nh}{\mu},0,0,0\right).$ (6)

####  2.3.2 Existence of Endemic Equilibrium (EE)

Using the definition of an equilibrium point and doing the substitutions the
endemic equilibrium, $E_{1}$, is expressed as

$E_{1}=(S^{*},E^{*},I^{*},R^{*}).$ (7)

where

$\left\\{\begin{array}[]{lll}S^{*}&=&\frac{(\mu+\varepsilon)(\mu+\tau+\gamma)}{\beta\varepsilon},\\\
&&\\\
E^{*}&=&\frac{(\mu+\rho)(\mu+\tau+\gamma)[hN\beta\varepsilon-\mu(\mu+\varepsilon)(\mu+\tau+\gamma)]}{\beta\varepsilon[(\mu+\rho)(\mu+\varepsilon)(\mu+\tau+\gamma)-\gamma\varepsilon\rho]},\\\
&&\\\
I^{*}&=&\frac{(\mu+\rho)[hN\beta\varepsilon-\mu(\mu+\varepsilon)(\mu+\tau+\gamma)]}{\beta[(\mu+\rho)(\mu+\varepsilon)(\mu+\tau+\gamma)-\gamma\varepsilon\rho]},\\\
&&\\\
R^{*}&=&\frac{\gamma[hN\varepsilon\beta-\mu(\mu+\varepsilon)(\mu+\tau+\gamma)]}{\beta[(\mu+\rho)(\mu+\varepsilon)(\mu+\tau+\gamma)-\gamma\varepsilon\rho]}.\end{array}\right.$
(8)

Therefore, the endemic equilibrium (EE), denoted $E_{1}$, given by equation
(7) exists whenever the associated reproduction threshold quantity, $R_{0}$,
exceeds unity. In other hand, no endemic equilibrium.

###  2.4 Basic Reproduction Number

The reproduction number, $R_{0}$, which measures the spread of infections in a
population, was computed by next generation matrix approach[13, 14].
Mathematically, $R_{0}$, is a treshold for stability of a disease-free
equilibrium (DFE) and is related to the peak and final size of an epidemic. In
other hand, $R_{0}$ is defined as the average number of new cases an infection
caused by one typical infected individual, in a population consisting of
susceptibles only [14].

Using this approach, the basic reproduction number, $R_{0}$, is given by
$R_{0}=\rho(-FV^{-1})$, where

$\begin{array}[]{ccc}F&=&\begin{pmatrix}0&\frac{\beta hN}{\mu}\\\ &\\\
0&0\end{pmatrix}\par\end{array}$ (9)

and

$\begin{array}[]{ccc}V^{-1}&=&\begin{pmatrix}\frac{-1}{\mu+\varepsilon}&0\\\
&\\\
\frac{-\varepsilon}{(\mu+\varepsilon)(\gamma+\mu+\tau)}&\frac{-1}{\gamma+\mu+\tau}\end{pmatrix}\end{array}.$
(10)

Which means that $R_{0}$ is given by dominant eigenvalues of $-FV^{-1}$, where
$F$ is the transmission part, describing the production of new infections, and
the $V^{-1}$ is the inverted matrix of $V$ where $V$ is the transition part,
describing changes in state. Thus,

$R_{0}=\frac{\beta Nh\varepsilon}{\mu(\mu+\varepsilon)(\gamma+\mu+\tau)}.$
(11)

###  2.5 Stability Analysis of disease equilibria

####  2.5.1 Stability of DFE

######  Theorem 2.1.

The disease free equilibrium $E_{0}$ of the model 1, given by equation 6, is
locally asymptotically stable (LAS) if $R_{0}\leq 1$ and unstable if
$R_{0}>1$.

Proof: This theorem will be proved based on the notions of the Jacobian
matrix. Let then be a point $M=(S,E,I,R)$ and $J(M)$ its Jacobian. We have:

$J(M)=\begin{pmatrix}-(\mu+\beta I)&0&-\beta S&\rho\\\ &&&\\\ \beta
I&-(\mu+\varepsilon)&\beta S&0\\\ &&&\\\ 0&\varepsilon&-(\mu+\tau+\gamma)&0\\\
&&&\\\ 0&0&\gamma&-(\mu+\rho)\end{pmatrix}.$ (12)

The Jacobian of DFE, $E_{0}$, expressed in equation 6. Finding the
eigenvalues, $z$, of $J(E_{0})$ amounts to solving the equation

$P(z)=0.$ (13)

where $P(z)$ is the characteristic polynomial of $J(E_{0})$, that is
$P(z)=|J(E_{0})-z\mathbb{I}_{4}|$.

The equation 13 becomes:

$(\mu+z)=0,$ (14)

or

$(\mu+\rho+z)=0$ (15)

or

$z^{2}+(2\mu+\varepsilon+\tau+\gamma)z+(\mu+\varepsilon)-\frac{h\beta\varepsilon
N}{\mu}=0$ (16)

Then, we have:

$z_{1}=-\mu,$ (17)

or

$z_{2}=-(\mu+\rho),$ (18)

or

$z^{2}+(2\mu+\varepsilon+\tau+\gamma)z+(\mu+\varepsilon)(\mu+\tau+\gamma)-\frac{h\beta\varepsilon
N}{\mu}=0,$ (19)

with Descartes’ rule of signs[15], if $R_{0}\leq 1$, all the coefficients of
the polyom characterizing the left side of the equation 19 are strictly
positive, then it does not have a positive root. From 17 18 19, if $R_{0}\leq
1$, $J(E_{0})$ has all its eigenvalues with strictly negative real part, hence
$E_{0}$ is locally asymptotically stable (LAS) by the Poincarré-Lyapunov
theorem[15, 16].

Otherwise, with the same rule, if $R_{0}>1$, there is a variation of
coefficients of $P(z)$, then $P(z)$ has at least one positive root. So,
$J(E_{0})$ has at least one eigenvalue with a strictly positive real part,
hence then $E_{0}$ is unstable by the Poincarré-Lyapunov theorem. $\Box$

The epidemiological implication of theorem 2.1 is that Tuberculosis spread can
be effectively controlled in the community when $R_{0}\leq 1$ that means if
the initial sizes of the populations of the model are in the basin of
attraction of the disease free equilibrium $E_{0}$. To ensure that elimination
of TB is independent of the initial sizes of the populations, it is necessary
to show that the DFE is globally asymptotically stable[17]. It is shown below
that tuberculosis will be eliminated from the community if the epidemiological
threshold can be reduced to a value below unity.

######  Theorem 2.2.

The disease free equilibrium $E_{0}$ of the model 1, given by equation 6, is
globally asymptotically stable (GAS) in $\Gamma$ whenever $R_{0}\leq 1$ and
unstable if $R_{0}>1$.

Proof: Consider the following Lyapunov candidate function for $E_{0}$:

$L=\varepsilon
E+(\mu+\varepsilon)I+AR~{}~{}\text{with}~{}A=\frac{\varepsilon\beta
N(h-\mu)}{\gamma\mu}.$ (20)

The first derivative of $L$ along the solutions of model 1 is

$\begin{array}[]{cll}\dot{L}&=&\varepsilon[\beta
SI-(\mu+\varepsilon)E]+(\mu+\varepsilon)[\varepsilon
E-(\mu+\tau+\gamma)I]+A[\gamma I-(\mu+\rho)R]\\\ &&\\\ &\leq&[\varepsilon\beta
S-(\mu+\varepsilon)(\mu+\tau+\gamma)+A\gamma]I\\\ &&\\\
&\leq&(\mu+\varepsilon)(\mu+\tau+\gamma)(R_{0}-1)I.\end{array}$ (21)

Since all the model parameters are non negative, it follows that $\dot{L}\leq
0$ for $R_{0}\leq 1$. For $R_{0}=1$, $\dot{L}=0$ if and only if $E=I=R=0$. As
$t\longrightarrow\infty$, substituting these values in model 1, we have that
$S\longrightarrow\frac{Nh}{\mu},~{}E\longrightarrow 0,~{}I\longrightarrow
0,~{}\text{and}~{}R\longrightarrow 0$. Hence, $L$ is a Lyapunov candidate
function on $\Gamma$, and the largest compact invariant set in
$\\{(S,E,I,R)\in\Gamma:\dot{L}=0\\}$ is the singleton $\\{E_{0}\\}$.
Therefore, by LaSalle’s Invariance Principle [17, 18], every solution of model
1, with initial conditions in $\Gamma$, approaches $E_{0}$ as
$t\longrightarrow\infty$, whenever $R_{0}\leq 1$ that is the disease free
equilibrium $E_{0}$ of the model 1, given by equation 6, is globally
asymptotically stable (GAS) in $\Gamma$ whenever $R_{0}\leq 1$ and unstable if
$R_{0}>1$. $\Box$

####  2.5.2 Stability of EE

######  Theorem 2.3.

The endemic equilibrium $E_{1}$ of the model 1, given by equation 7, is
locally asymptotically stable (LAS) if $R_{0}>1$.

Proof: From $N=S+E+I+R$, we deduce that $R=N-S-E-I$. By reducing dimension of
system 1 while keeping the variables $S,E,I$, we have:

$\left\\{\begin{array}[]{ccc}\frac{dS(t)}{dt}&=&hN(t)-\mu S(t)-\beta
S(t)I(t)+\rho(N-S-E-I),\\\ &&\\\ \frac{dE(t)}{dt}&=&\beta
S(t)I(t)-(\mu+\varepsilon)E(t),\\\ &&\\\ \frac{dI(t)}{dt}&=&\varepsilon
E(t)-(\mu+\tau+\gamma)I(t).\end{array}\right..$ (22)

Let $W=(S,E,I)$ be a point and $J(W)$ its Jacobian. We have:

$J(W)=\begin{pmatrix}-(\mu+\rho+\beta I)&-\rho&-(\rho+\beta S)\\\ &&&\\\ \beta
I&-(\mu+\varepsilon)&\beta S\\\ &&&\\\ 0&\varepsilon&-(\mu+\tau+\gamma)\\\
\end{pmatrix}.$ (23)

In this form, the reduction of endemic equilibrium point EE, $E_{1}$, becomes
$E_{1}=(S^{*},E^{*},I^{*})$. Let us look for the eigenvalues of $J(E_{1})$,
this amounts to find the roots of the characteristic polynomial $P(z)$ that
means resolve $P(z)=0$. By factoring determinant above, we have

$P(z)=A(z)B(z)C(z)$ (24)

with

$A(z)=-[(\mu+\rho+\beta I^{*})+z],$ (25)

$B(z)=z^{2}+(2\mu+\varepsilon+\rho+\beta I^{*})z+(\mu+\rho+\beta
I^{*})(\mu+\varepsilon)+\rho\beta I^{*},$ (26)

and

$\begin{array}[]{lll}C(z)&=&z^{4}+(3\mu+2\beta
I^{*}+2\rho+\varepsilon+\tau+\gamma)z^{3}\\\ &&\\\ &&+[(2\mu+\rho+\beta
I^{*}+\tau+\gamma)(2\mu+\varepsilon+\rho+\beta I^{*})\\\ &&\\\
&&+(2\mu+\rho+\beta I^{*})(\mu+\varepsilon)+(\mu+\rho+\beta
I^{*})(\mu+\tau+\gamma)+\rho\beta I^{*}]z^{2}\\\ &&\\\ &&+[(\mu+\rho+\beta
I^{*})(\mu+\tau+\gamma)(2\mu+\varepsilon+\rho+\beta I^{*})+\varepsilon\beta
S^{*}]z\\\ &&\\\ &&+[(\mu+\rho+\beta I^{*})(\mu+\varepsilon)+\rho+\beta
I^{*}][2\mu+\rho+\beta I^{*}\tau+\gamma\\\ &&\\\ &&+(\mu+\rho+\beta
I^{*})(\mu+\rho+\beta I^{*})(\mu+\tau+\gamma)]\\\ &&+\varepsilon\beta
S^{*}(\mu+\rho+\beta I^{*})+\rho\beta I^{*}(\rho+\beta S^{*})\end{array}.$
(27)

Resolve $P(z)=0\Longrightarrow
A(z)=0,~{}~{}\text{or}~{}~{}B(z)=0,~{}~{}\text{or}~{}~{}C(z)=0$. Thus, with
Descartes’ rule of signs[15], if $R_{0}>1$, all the coefficients of the polyom
characterizing the right side of the equations 26, 27 are strictly positive,
then they do not have a positive root. From 25 26 27, if $R_{0}>1$, $J(E_{1})$
has all its eigenvalues with strictly negative real part, hence $E_{1}$ is
locally asymptotically stable (LAS) by the Poincarré-Lyapunov theorem[15, 16].
Therefore, this is the end of the theorem 2.3. $\Box$

This implies that the tuberculosis persist and invade the population if
$R_{0}>1$.

##  3 Numerical simulation

In this section, we are going to demonstrate the behavior of the proposed TB
infection model through numerical simulations while comparing the model with
actual tuberculosis data in Burundi. Data from Burundi were collected on
quarterly (resp. annually) from the Ministry of Public Health between 2011 and
2019 (resp. between 1997 and 2020). As shown in Table 1, the parameter values
for simulation are either from literature or reasonably chosen estimates. We
used Python to run simulations and the Runge-Kutta method was used to
approximate solutions of the SEIR-system.

Table 1: Parameters values Parameter |  Parameter Value and references
---|---
h |  Human birth rate is 38.377 per year per 1000 inhabitants [16]
$\beta$ |  0.9 Probability of being infected for a susceptible meted by an infectious [19]
$\rho$ |  0.897 Rate of progression from recovered to susceptible state Estimated
$\gamma$ |  0.525 Rate of progression from infectious to recovered state Estimated
$\mu$ |  Human (natural) mortality rate is 7.766 per year per 1000 inhabitants [16]
$\varepsilon$ |  0.25 Rate of progression from exposed to infectious state [20]
$\tau$ |  21.0 (per 100,000 people) in 2019 human tuberculosis-induced death rate [21]

###  3.1 Fitting the model to tuberculosis data in Burundi

In this paragraph, we use systems 1 and parameters listed in the table (1) to
fit the tuberculosis model to the annual data of population which are
diagnostics of tuberculosis in Burundi from $1985$ to $2020$. So,
$S_{0}=8,895,836,~{}E_{0}=375,~{}I_{0}=1789$ and $R_{0}=60000$ are assumed
initial conditions of states

Figure 2: Fitting the model to tuberculosis data in Burundi

The yearly data in Figure (2) follows the trend more closely. In this case, we
have obtained a good fit. Additionally, the yearly increase in TB cases could
be attributed to poor interventions implementation.

###  3.2 Simulation of TB transmission dynamics

Figure (3) demonstrates that if $R_{0}$ is more than 1, each infectious person
spreads the disease to more than one new person throughout the time of
contagion, which results in the spread of the disease throughout the
population and the emergence of an epidemic

Figure 3: Human population dynamics in relation to infectious, and recovered
humans when $R_{0}>1$

and Figure (4) shows that if $R_{0}$ less than 1 suggests that each person
infects, on average, less than one new person, resulting in disease
extinction.

Figure 4: Human population dynamics in relation to infectious, and recovered
humans when $R_{0}<1$

##  4 Sensitivity Analysis

In order to determine the relative importance of model parameters on the
disease infection, a sensitivity of parameters of the model system 1 is
carried out. Sensitivity indices allow us to measure the relative change in a
state variable when a parameter changes [22]. The numerical calculation of the
sensitivity indices also allows the determination of parameters which have a
reasonable impact on the basic reproduction number $R_{0}$ and which of the
parameters is most sensitive which can help in the eradication of the disease
in the population [23].

The normalized forward sensitivity index of a variable to a parameter is the
ratio of the relative change in the variable to the relative change in the
parameter. When the variable is a differentiable function of the parameter,
the sensitivity index may be alternatively defined using partial derivatives
[24, 22, 23].

######  Definition 4.1.

The normalized forward sensitivity index of a variable, $u(p)$, that depends
differentiably on a parameter, $p$, is defined as:

$K_{p}^{u}=\frac{\partial u}{\partial
p}\times\frac{p}{u},~{}~{}~{}\text{with}~{}~{}u\neq 0$ (28)

Therefore, we derive an analytical expression for the sensitivity of $R_{0}$
with $R_{0}$ given by the expression 11. Thus,

$K_{p_{i}}^{R_{0}}=\frac{\partial R_{0}}{\partial
p_{i}}\times\frac{p_{i}}{R_{0}},~{}~{}~{}\text{with}~{}~{}p_{i}:i=1,\cdots,n$
(29)

$p_{i}$ denotes each parameter involved in $R_{0}$.

According to Table 1, and the equation 29, we plot sensitivity index of each
parameter with respect to the $R_{0}$ given by equation 11. Sensitivity
indices on $R_{0}$ are showed on the the following bar chart

Figure 5: Bar chart for sensitivity

All parameters of the model are assumed to be non-negative, let the total
population be 100,000. Hence most sensitive parameter is being the highest
positive index. On Figure 5 we can list
$\beta,~{}h~{}\text{and}~{}\varepsilon$. This implies that increasing these
parameters causes $R_{0}$ to increase by $100\%$. Thus, since $R_{0}$
continues to be higher, the epidemic of disease infection tends to occur. On
the other hand, the parameters whose sensitivity is negative
($\gamma,~{}\tau~{}\text{and}~{}\mu$) have a highest negative impact on
$R_{0}$. For the following, a descriptive analysis as well as the linear
adjustment of our model will be carried out.

##  5 Conclusion

In this paper, a deterministic mathematical model for the transmission
dynamics of TB in Burundi has been carried out to study the stability of both
disease-free and endemic equilibrium point. Using matrix generation approach,
the basic reproduction number $R_{0}$ was computed. Therefore, the disease-
free equilibrium of the model obtained is both locally and globally stable for
$R_{0}\leq 1$. It is also shown that the endemic equilibrium solution of the
model is globally asymptotically stable if $R_{0}>1$. Sensitivity Analysis has
showed that $\beta,~{}h~{}\text{and}~{}\varepsilon$ are being the highest
positive index.

Using data from Ministry of Public Health and the Fight against AIDS,
parameters of the model have been identified and the model is shown to
describe TB dynamic in Burundi from 1997 to 2019 (annual data) and from 2011
to 2019 (quarterly data).

The results of this research support the idea that in sub-Saharan Africa
people in general and in Burundi in particularly should be strongly encouraged
to seek a diagnosis of TB, and the success rate the treatment must be
correlated with the population of infectious diagnosed. Then, the population
of diagnosed infectious diseases and the number of TB-related deaths will
decrease.

## References

  * [1] TM Daniel.  Historical review.  The history of tuberculosis, in respiratory medicine, 100(11):1862–1870, 2006.
  * [2] Alimuddin Zumla, Andrew George, Virendra Sharma, Rt Hon Nick Herbert, Aaron Oxley, and Matt Oliver.  The who 2014 global tuberculosis report—further to go.  The Lancet Global Health, 3(1):e10–e12, 2015.
  * [3] A Mandal.  Diarrhea epidemiology (online) diakses dari http://www. newsmedical. net/health.  Diarrhea-Epidemiology. aspx pada, 18, 2013.
  * [4] Marilyn Ronoh, Rym Jaroudi, Patrick Fotso, Victor Kamdoum, Nancy Matendechere, Josephine Wairimu, Rose Auma, and Jonnes Lugoye.  A mathematical model of tuberculosis with drug resistance effects.  Applied Mathematics, 7(12):1303–1316, 2016.
  * [5] Graham AW Rook and Barry R Bloom.  Mechanisms of pathogenesis in tuberculosis.  Tuberculosis: pathogenesis, protection, and control, pages 485–501, 1994.
  * [6] Bess Miller and Dixie E Snider Jr.  Physician noncompliance with tuberculosis preventive measures, 1987.
  * [7] Carlos Castillo-Chavez, Zhilan Feng, et al.  Mathematical models for the disease dynamics of tuberculosis.  1996\.
  * [8] Wikipedia.  https://en.wikipedia.org/wiki/Isoniazed, 28th september, 2021.
  * [9] K6572870412 Dietz.  Epidemiologic interference of virus populations.  Journal of mathematical biology, 8:291–300, 1979.
  * [10] KS Warren, RM Anderson, V Capasso, AD Cliff, K Dietz, F Fenner, RN TW-Fiennes, Z Grossman, H Knolle, PG Mann, et al.  Population biology of infectious diseases, 1982.
  * [11] Zhilan Feng and Horst R Thieme.  Recurrent outbreaks of childhood diseases revisited: the impact of isolation.  Mathematical biosciences, 128(1-2):93–130, 1995.
  * [12] Hans Waaler, Anton Geser, and Stig Andersen.  The use of mathematical models in the study of the epidemiology of tuberculosis.  American Journal of Public Health and the Nations Health, 52(6):1002–1013, 1962.
  * [13] Hyun Mo Yang.  The basic reproduction number obtained from jacobian and next generation matrices–a case study of dengue transmission modelling.  Biosystems, 126:52–75, 2014.
  * [14] Pauline Van den Driessche and James Watmough.  Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission.  Mathematical biosciences, 180(1-2):29–48, 2002.
  * [15] Mamadou Lamine Diouf.  Analyse de modèles épidémiologiques à plusieurs classes d’infectés: stabilité et observabilité.  PhD thesis, Université Gaston BERGER de Saint-Louis du Sénégal, 2016.
  * [16] Egide Ndamuzi and Paterne Gahungu.  Mathematical modeling of malaria transmission dynamics: Case of burundi.  Journal of Applied Mathematics and Physics, 9(10):2447–2460, 2021\.
  * [17] Junli Liu and Tailei Zhang.  Global stability for a tuberculosis model.  Mathematical and Computer Modelling, 54(1-2):836–845, 2011.
  * [18] Joseph P La Salle.  The stability of dynamical systems.  SIAM, 1976.
  * [19] Jannatun Nayeem, Israt Sultana, et al.  Mathematical analysis of the transmission dynamics of tuberculosis.  American Journal of Computational Mathematics, 9(03):158, 2019.
  * [20] Ihsan Ullah, Saeed Ahmad, Qasem Al-Mdallal, Zareen A Khan, Hasib Khan, and Aziz Khan.  Stability analysis of a dynamical model of tuberculosis with incomplete treatment.  Advances in Difference Equations, 2020(1):1–14, 2020.
  * [21] atlas.  https://knoema.fr/atlas/Burundi, 30th october, 2021.
  * [22] Nakul Chitnis, James M Hyman, and Jim M Cushing.  Determining important parameters in the spread of malaria through the sensitivity analysis of a mathematical model.  Bulletin of mathematical biology, 70:1272–1296, 2008.
  * [23] TS Faniran, AO Falade, and J Ogunsanwo.  Sensitivity analysis of an untreated liver-stage malaria.  Journal of Computer Science & Computational Mathematics, 10:61–68, 2020.
  * [24] Titus Orwa, Rachel Mbogo, and Livingstone Luboobi.  Mathematical model for the in-host malaria dynamics subject to malaria vaccines.  Letters in Biomathematics, 5(1):222–251, 2018.


[1, $\star$]AngelinaAgabin [2,3 $\star$]J. XavierProchaska

1]Department of Applied Math, University of California, Santa Cruz, CA, 95064,
USA 2]Affiliate of the Department of Ocean Sciences, University of California,
Santa Cruz, CA, 95064, USA 3]Department of Astronomy and Astrophysics,
University of California, Santa Cruz, CA, 95064, USA

# Reconstructing Sea Surface Temperature Images: A Masked Autoencoder Approach
for Cloud Masking and Reconstruction

###### Abstract

This thesis presents a new algorithm to mitigate cloud masking in the analysis
of sea surface temperature (SST) data generated by remote sensing
technologies, e.g., the Level-2 Visible-Infrared Imager-Radiometer Suite
(VIIRS). Clouds interfere with the analysis of all remote sensing data using
wavelengths shorter than $\approx 12$ microns, significantly limiting the
quantity of usable data and creating a biased geographical distribution
(towards equatorial and coastal regions). To address this issue, we propose an
unsupervised machine learning algorithm called Enki which uses a Vision
Transformer with Masked Autoencoding to reconstruct masked pixels. We train
four different models of Enki with varying mask ratios (referred to as $t$) of
10%, 35%, 50%, and 75% on the generated Ocean General Circulation Model (OGCM)
dataset referred to as LLC4320. To evaluate performance, we reconstruct a
validation set of LLC4320 SST images with random “clouds” corrupting $p$=10%,
20%, 30%, 40%, 50% of the images with individual patches of $4\times 4\,\rm
pixel^{2}$. We consistently find that at all levels of $p$ there is one or
multiple models that reconstruct the images with a mean RMSE of less than
$\approx 0.03$K, i.e. lower than the estimated sensor error of VIIRS data.
Similarly, at the individual patch level, the reconstructions have RMSE
$\approx 8\times$ smaller than the fluctuations in the patch. And, as
anticipated, reconstruction errors are larger for images with a higher degree
of complexity. Our analysis also reveals that patches along the image border
have systematically higher reconstruction error; we recommend ignoring these
in production. We conclude that Enki shows great promise to surpass in-
painting as a means of reconstructing cloud masking. Future research will
develop Enki to reconstruct real-world data.

## Plain Language Summary

This thesis presents a solution to the problem of cloud masking in the
analysis of sea surface temperature (SST) data generated by remote sensing
technology. With the vast amount of SST data generated each year, unsupervised
machine learning algorithms have proven effective in identifying patterns
within the dataset. However, cloud coverage can interfere with the analysis of
the data as convolutional neural network models like Ulmo are not equipped to
handle masked images. To address this issue, we propose an unsupervised
machine learning algorithm, called Enki which uses a Vision Transformer with
Masked Autoencoding to reconstruct pixels that are masked out by clouds.

##  1 Introduction

Remote sensing technology has revolutionized our ability to monitor and
understand Earth’s systems. Sensors like the Visible Infrared Imaging
Radiometer Suite (VIIRS; Jonasson and Ignatov, 2019) or the Moderate
Resolution Imaging Spectroradiometer (MODIS;
https://oceancolor.gsfc.nasa.gov/data/aqua/) are used for tasks such as
weather forecasting, monitoring climate change, atmospheric studies and more.
One of the most important oceanographic variables to observe is sea surface
temperature (SST), which impacts and or tracks key aspects of our ecosystem,
e.g., marine dynamics, weather patterns, global warming (Isern-Fontanet et
al., 2017). VIIRS and MODIS both capture SST data at high-spatial resolution
($\approx 1$ km) with global coverage and twice-daily measurements for the
past one (VIIRS) or two (MODIS) decades.

Sensors like VIIRS have enabled routine monitoring of the ocean, but because
remote sensing technology generates such large datasets (e.g., the National
Oceanic and Atmospheric Administration’s Level-2P (L2P), 2nd full-mission
reanalysis (RAN2) of the VIIRS dataset is nearly 100Tb in size), there is a
growing need for automated and efficient ways to sort through and analyze the
data. Unsupervised machine learning (ML) is one approach that has shown great
promise in this regard. Specifically, by applying unsupervised machine
learning to remote sensing data, we can extract meaningful information from
large datasets without the need for human intervention or (expensive)
labeling. This enables one, for example, to observe ocean patterns that may
reveal the impact of global warming or to identify shifting currents,
upwelling or marine heat waves that impact human activity (Prochaska et al.,
2023; Levy et al., 2018; Penven et al., 2006; Prochaska et al., 2023). Within
Physical Oceanography, such models may identify patterns or outliers in sea
surface temperature data, such as the development of oceanic fronts or eddies
(McWilliams, 2017; Chelton et al., 2011).

With the rise of ML, a type of artificial intelligence (AI), analysis of
remote sensing data, a difficult challenge inherent to these data rises to the
fore: corrupt data. Clouds, the most impactful corruption, stymie remote
sensing observations leading to masking and missing data on nearly all spatial
scales. These are identified with custom algorithms and flagged in the Level 2
products produced by NASA and NOAA among others (Wu et al., 2017). Corrupt
data can greatly reduce the application and analysis of remote sensing data
(Figure 1), especially in select geographical regions (Figures 2). The problem
is especially acute for many standard deep-learning AI models (e.g.
convolutional neural networks, CNNs) which generally cannot do not allow for
missing or masked data.

Figure 1: Cumulative number of cutouts $N_{c}$ that are available for analysis
as a function of cloud coverage $\rm CC$. At small values (${\rm CC}<10\%$),
there is a greater than exponential rise in the amount of data available for
analysis. Increasing from ${\rm CC}=5\%$ to 20% would make $\approx 3\times$
more data available and one would gain an order-of-magnitude by extending to
${\rm CC}=40\%$.  Figure 2:  Geographical distribution of the data with $\leq
5\%$ masked data (primarily clouds) in the nighttime L2 MODIS dataset of SST
imagery (see Prochaska et al., 2021, for full details). Clearly, limiting to
nearly cloud-free data greatly limits SST analysis in various regions. Note
the color bar is on a log10 scale.

For example, consider the unsupervised ML model Ulmo by our team (Prochaska et
al., 2021), a probabilistic autoencoder architecture originally trained on SST
images that span 150x150km2 regions, which this paper will refer to as
cutouts, from the L2 SST Aqua MODIS dataset, and later retrained on RAN2 VIIRS
SST data (Gallmeier et al., 2023). Ulmo utilizes two distinct algorithms: an
autoencoder and a normalizing flow. The autoencoder reduces the dimensionality
of the original image into a latent vector, which represents the model’s
understanding of the image. These latent vectors are then processed through a
normalizing flow to convert them into a single variable known as the log
likelihood ($\rm LL$). The lower the LL, the more uncommon the image and the
greater its complexity Prochaska et al. (2021). This allows one to perform
outlier detection to find interesting patterns.

However one major limitation of Ulmo and CNNs in general is that the input
images need to be “whole”. For SST measurements derived from near-infrared
wavelengths, therefore, clouds pose a major problem. In the case of Ulmo, the
authors primarily trained on cutouts that are cloud-free or have low cloud
coverage (${\rm CC}\leq 5\%$), and used a standard inpainting algorithm to
fill in the masked pixels. This resulted in a biased analysis of the ocean
with the probability of a cutout decreasing away from the continental
boundaries except in the Equatorial region (Figure 2), and a significantly
reduced fraction of the total data available (Figure 1).

Figure 1 highlights the potential gains of cloud mititgation. Even extending
to data from ${\rm CC}=5\%$ to ${\rm CC}=20\%$ would increase the dataset by
$\sim 3\times$ while extending to 40% would increase it by over one order of
magnitude. However, our experiments with the inpainting algorithm used in the
Ulmo (Navier-Stokes Bertalmio et al., 2001) and subsequent analyses indicate
poor performance at ${\rm CC}>5\%$, and even reveal important systematic
errors at lower cloud coverage fractions.

To make progress on this critical issue, we propose a new solution – use a
novel, ML algorithm based on natural language processing (NLP) to reconstruct
masked data. Our model, named Enki, is a Masked Autoencoder built upon the
architecture of a Vision Transformer (ViTMAE) trained on LLC SST cutouts (He
et al., 2021). By reconstructing the masked out pixels, this allows for SST
cutouts to be analyzed by algorithms like Ulmo, and more broadly, provide
accurate estimates for missing data for other scientific and commercial
applications.

This Masters thesis to the Department of Scientific Computing and Applied
Maths at the University of California, Santa Cruz describes the development,
testing, and first results from Enki. The thesis is organized as follows:
Section 2 covers the architecture of ViTMAE. Section 3 goes over the
environment in which Enki was trained, the specifics of our training dataset
and validation set, and the hyperparameters we used during training. Section 4
goes over the qualitative and quantitative results as we examine Enki ’s
reconstructions.

Figure 3: Architecture of our ViTMAE named Enki. In this example, the image is
64x64 pixels and broken up into patch sizes of 16x16 pixels with 40% of the
image masked out. An image is broken down into patches and masked. The
unmasked patches are flattened and embedded by a linear projection with
positional embeddings which are then run through the encoder, returning the
encoded patches which are then run through the decoder along with the masked
(unfilled) tokens. This returns another latent vector, and another linear
projection layer outputs this vector as an image with the same dimensions of
the original image. The model of Enki we train reconstructs cutouts of size
$64\times 64$ pixel2 with a patch size of $4\times 4$ pixel2 and a latent
vector size of 256.

##  2 Architecture

Enki is a Masked Autoencoder (MAE) neural network with a Vision Transformer
(ViT) architecture (He et al., 2021). As an autoencoder, it has two primary
parts: (1) the encoder and (2) the decoder (e.g. Baldi, 2012). The encoder
specializes in transforming an input image into a lower-dimensional latent
space, known as a latent vector. In machine learning, latent vectors are
reduced, numerical representations of data which (ideally) capture the
essential characteristics or features of the data. The decoder transforms this
latent vector back into a reconstructed output image with the same dimensions
of the input space.

Transformers are another type of neural network that also have an encoder and
decoder type architecture but differ from autoencoders in their structure and
use. Unlike autoencoders, transformers tokenize inputs. Each token contains a
latent representation, and the transformer applies self-attention to weigh the
contribution of each latent vector to the final data representation. Self-
attention works by computing a set of attention weights for each latent
vector, based on its similarity and association to other latent vectors in the
image. Transformers are primarily used in NLP for tasks such as language
translation. An example would be if a transformer was given an English
sentence with the goal of reconstructing the sentence in Japanese. The
transformer tokenizes each word, performs self-attention to find the relations
between the words and produces the latent vector, and then reconstructs the
sentence in Japanese from the latent vector.

A ViT is a Transformer encoder applied to images (Dosovitskiy et al., 2020).
Similar to a sentence, the image is tokenized by breaking it into non-
overlapping patches, each of which is assigned a positional embedding. These
are then fed through a standard Transformer encoder. Self-attention is also
performed on the patches to weigh contributions of each patch to the final
image representation. ViT is primarily used for image classification (i.e.,
without a Transfomer decoder) and is therefore a supervised learning model
requiring labeled data for training.

He et al. (2021) developed a masked autoencoder based on the ViT which
utilizes both the encoder and decoder and learns to reconstruct images with
masked/missing data. The masked nature allows for unsupervised training which
is ideal for handling unlabeled big data. ViTMAE works by taking an image,
breaking it up into patches, masking out a portion of the patches, and then
tokenizing them. These patches are then fed through a transformer (the
encoder), and a latent vector is returned. The latent vector is then given to
another transformer block (the decoder), which reconstructs the image and a
linear projection layer converts it back into the original image’s dimension
size.

Our architecture, named Enki and described in Figure 3, follows closely that
of the original ViTMAE. We changed the hyper parameters of image size to fit
our $64\times 64$ pixel2 cutouts, and reduced the patch size down to $4\times
4$ pixel2. A small patch size is preferred because real clouds are not large
squares. A future version of Enki will convert cloud mask into a set of
patches, and we can reproduce a mask with squares while minimizing the number
of good pixels masked.

On the other hand, smaller patch sizes are more computationally expensive.
Decreasing the patch size down to even $2\times 2$ pixel2 will increase the
total number of patches from 256 to 1024 with an increase of $4\times$ in
compute.

Config | Value
---|---
Models and Masking Ratios | t=10, 35, 50, 75
Base Learning rate (t=10, 35, 50) | 1e-4
Base Learning rate (t=75) | 1.5e-4
Weight Decay | 0.05
Warmup epochs | 40
Training Epochs | 400
Batch size | 64
Image size | (64, 64, 1)
Patch size | (4, 4)
Encoder Dimension | 258
Decoder Dimension | 512
Table 1: The training hyperparameters for the model are summarized in this
table, highlighting the different learning rates chosen for each specific
model. Additionally, a 256 embedding dimension was initially selected with the
intention of utilizing the latent vectors for training Ulmo. However, it was
discovered that this approach was not feasible due to the computationally
expensive nature of training Ulmo with each patch generating a 256-sized
latent vector. As a result, alternative strategies were explored to address
this limitation. Figure 4: A visual example of the reconstruction process. 1)
We start with the Original image; 2) We convert this image into patches, and
mask out a $p$ percentage of the image. In this example, we use $p=75$; 3) The
unmasked patches are run through a transformer block (the encoder; see Figure
3) that returns latent representations of the patch. Each unmasked patch has a
latent vector of size 256; 4) The full set of patches (latent vectors of
unmasked patches + empty masked patches) are now run through another
transformer block (the decoder) that reconstructs the missing patches. This
image is converted back into the original dimension space of the original
image. However, notice that what the decoder returns has notable patches where
the image was not masked. We suspect this is because when the model is
evaluating MSE, it only evaluates the reconstructed masked pixels. This means
when reconstructing, the model does not care what the unmasked patches look
like; and 5) We take the unmasked patches from the Original image and replace
the missing patches with the reconstructed patches from the decoder to create
our final reconstructed image.

##  3 Methodology

Enki is trained in two steps: pre-training and fine-tuning. Pre-training is
typically used to train a deep-learning model on a general corpus, while fine-
tuning is used to adapt the model to a specific dataset. A model of ViTMAE
could be trained on a general corpus like ImageNet-1K and then fine-tuned on
another dataset to specialize in reconstructing birds. However, when a model
of ViTMAE is trained, the hyperparameters used during training cannot be
changed when reconstructing new images outside of training. For example, a
model of ViTMAE was trained on 128x128, three color channel images with a
16x16 patch size, must be applied to data with the same shape. Our training
set did not match the parameters of the available trained models. Therefore
Enki was pre-trained from scratch to suit our needs. Pre-training for Enki was
performed using a modified version of Facebook AI Research group’s code to run
in Kubernetes with an h5 file containing our training dataset. Multiple models
of Enki were trained on different masking ratios ($t$) to study how its affect
on Enki performance for image reconstruction.

Finally, we have experimented with different masking ratios during and after
pre-training. We distinguish these two values as the masking ratio a model was
trained at as $t$ (e.g. a model trained at 10% masking is referred to as
$t$=10) and the masking ratio for the input image for reconstruction is
referred to as $p$ (e.g. an image being reconstructed at 10% masking is
referred to as $p$=10).

###  3.1 Training Dataset

In unsupervised learning, the dataset does not have labels. Our training
dataset is comprised of millions of SST images for Enki to learn patterns. For
this project, the optimal dataset would be one that is free of clouds and
noise because we do not want Enki to learn to reconstruct unmasked clouds or
add noise into reconstructed patches(Wu et al., 2017). The dataset should
adequately capture both typical and uncommon SST dynamics with a uniform
coverage of the ocean.

For pretraining, we used a dataset generated from a Ocean General Circulation
Model (OGCM) from the Estimating the Circulation and Climate of the Ocean
(ECCO) project in a collaborative effort between the Massachusetts Institute
of Technology (MIT), the Jet Propulsion Laboratory (JPL), and the NASA Ames
Research Center (ARC). Specifically, we use the $\frac{1}{48}^{\circ}$,
90-level simulation known as LLC4320. The LLC4320 model outputs do not contain
clouds (only a statistical model of the atmosphere) or sensor noise.

Furthermore, while we have large sets of cutouts of true data from MODIS and
VIIRS that are nominally “cloud free”, our experience is that many of these
contain clouds that were missed by standard screening algorithms (see
Gallmeier et al., 2023, for examples). Working with a cloud-free VIIRS dataset
may also bias Enki towards behaving best in cloud-free regions (e.g. Figure
2). LLC4320 data does not suffer from this issue; we can generate a training
set with uniform coverage of the ocean in both location and time.
Additionally, Gallmeier et al. (2023) shows that the distribution of sea
surface temperature anomaly (SSTa) patterns present in the VIIRS observations
are generally well-predicted by the LLC model. This provides confidence that a
ViTMAE trained on LLC4320 outputs may well reflect patterns within the real
ocean.

The training dataset for Enki consists of LLC4320 cutouts from 2011-11-17 to
2012-11-15, inclusive. LLC4320 cutouts were extracted to cover $144\times 144$
km2 sampling and resized using the local mean to $64\times 64$ pixel2 cutouts.
The geographical distribution uniformly samples the ocean between latitudes of
90S and 57N. Real satellite data contain noise from the atmosphere and from
instruments, but because the goal of reconstruction is not to reconstruct the
noise in images, we did not impute noise in the LLC4320 cutouts (in contrast
to Gallmeier et al., 2023). In total, the training dataset ends up being $\sim
2.6$ million cutouts. Although it is possible to generate larger datasets for
training purposes, we encountered resource limitations during the training
process, and expanding the dataset size would significantly extend the
training duration. Unlike datasets such as ImageNet, which contain a diverse
range of objects like animals, plants, and various objects, SST patterns
exhibit less variability. Given this characteristic, we hoped that the
inclusion of 2.6 million cutouts in our dataset adequately captures the
dynamics present in the ocean.

Figure 5: Training loss of the $t=75$ Enki model. The lost function
implemented was the MSE of the reconstructed masked out patches. As typical of
deep-learning models, there is rapid improvement in the loss at early stages
followed by a shallow plateau at late times.

###  3.2 Training

We performed a hyperparameter scan of the masking ratio during training $t$.
Four different models of Enki were trained with 10%, 35%, 50%, and 75%
masking. These variations are referred to as $t=10,35,50,75$ respectively.
Each model was trained for 400 epochs on the complete dataset of 2.6 million
LLC4320 SST cutouts. During training, several models errored if the learning
rate was too high. For these models, we adjusted the learning rates from the
default value; the $t=75$ model was trained at 0.00015, while the $t=10,35,50$
models were trained at 0.0001. During a single training epoch, Enki iterates
through every cutout in the dataset and masks out the specified masking ratio
with random $4\times 4\,\rm pixel^{2}$ patches. Performance was evaluated
using the default Mean Squared Error (MSE) loss function, but going forward we
may consider implementing other loss functions.

Every epoch of training uses new randomized masks for each cutout. This forces
the model to see significant variations and better generalize. All models were
trained in the Nautilus Pacific Rim Platform with Kubernetes using 4 CPUs,
75Gb of RAM, and 8 NVIDIA-A10 GPUs. These specs were chosen after profiling
training speeds with available resources. Training times varied based on
masking ratios. Higher masking ratios completed in shorter training times,
with the $t=75$ model training the fastest ($\approx 22$ minutes per epoch).
Lower masking ratios require more computations to create the latent vectors;
the $t=10$ has a 40% longer training time ($\approx 30$ minutes per epoch).

###  3.3 Validation Dataset

In unsupervised learning, the validation set is a subset of the data not used
during training. It is used to evaluate the performance of a model with data
the model has not yet seen. The validation set therefore assesses model
performance while avoiding any effects of overfitting during training.
Overfitting refers to when a model learns the dataset too well that it is able
to give accurate predictions, or in our case reconstructions, for the training
dataset, but has not generalized well to reconstruct new data. While 2.6
million cutouts is a substantial dataset, it is still small compared to the
100Tb of data in the VIIRS dataset. If we want Enki to perform well, we wish
it to reconstruct data outside of the training set.

The validation set we created for this project was a batch of 680,000 LLC4320
cutouts from the same locations of the training dataset but offset by $\approx
2$ weeks from the training set. These cutouts follow similar cloud free and no
noise and were drawn uniformly over the ocean. While not within the scope of
this thesis, the validation set for future studies could contain VIIRS or
MODIS cutouts to determine Enki’s effectiveness reconstructing real cutouts.

##  4 Analysis and Results

Following the completion of training, we used the validation set to
reconstruct cutouts with the four Enki models trained at $t=10,35,50,75$. For
each model, we reconstructed images at a masking ratio of $p=10,20,30,40,50$.
We chose to limit our range to $p\leq 50$ because further masking would run
the risk of Enki having insufficient information to interpolate accurately.
The seed for the masking was not set, so masks for all the reconstructions
differ. Reconstructed images were created by replacing masked patches with
reconstructed patches. The results are invariably imperfect and we anticipate
larger errors at higher masking ratios. To assess the magnitude of the error,
we calculate the Root Mean Squared Error (RMSE):

${RMSE}(y,\hat{y})=\sqrt{\frac{\sum_{i=0}^{N-1}(y_{i}-\hat{y}_{i})^{2}}{N}}$
(1)

where $y$ represents the original image, $\hat{y}$ represents the
reconstructed image, and $N$ represents the number of pixels that were masked
out. Throughout, the RMSE is only calculated on masked pixels $i$.

Analysis of the RMSE is performed at the patch level and on images as a whole
(excluding the edge patches, which we will discuss in depth in Section 4.2).
As discussed in Section 4.3, we identified a bias for select $t$,$p$ values.
In nearly all cases, we correct for this bias before evaluating the RMSE. In
the following sub-sections, we present results including the performance of
the Enki models with various $t$ values.

Figure 6:  Qualitative assessment of the reconstruction for a single image
with Enki. This high-complexity image (${\rm LL}\approx-19000$) was imputed
with a $p=50$ mask and we applied the $t=50$ Enki model. Despite the high
masking ratio, the reconstruction reproduces the original at high fidelity.
The primary exception are features along the sharpest gradient (residuals of
7.3 and -7.5) and the upper-right corner (residuals of -2.5) where the model
is forced to extrapolate.

###  4.1 Qualitative Results

Figure 7:  Galleries of reconstructions by the different models. In each case,
the cutout is reconstructed using the Enki-trained model that matches the
masking ratio $p$. A collection of reconstruction galleries showcasing the
results obtained from various models is presented in this study. Each gallery
exhibits the reconstruction of a different cutout being reconstructed by a
model at the mask training ratio they were trained at. Different LL were
intentionally selected to demonstrate how Enki performs across a spectrum of
ocean dynamics.

Consider first a visual inspection and assessment of the image reconstructions
by Enki. Figure 6 shows a representative example for a higher complexity SST
cutout. The example has a $p=50$ mask ratio and we implemented the $t=50$ Enki
model.

Overall, the Enki model demonstrates strong performance in capturing the
underlying structure of the SST patterns, even at higher masking ratios. The
reconstructed images exhibit smooth transitions between the unmasked and
reconstructed patches and no significant offsets, eliminating the need for
bias correction. However, some residual errors are observed along sharper
gradients, where the model accurately estimates the border between different
values but falls short of perfection. We also see that in cases where certain
features, such as the curved point of negative (blue) values, are completely
masked out, Enki produces sensible, albeit imperfect, reconstructions. A
similar occurrence can be observed in the upper-right corner where the
gradient is mostly masked out, and the reconstructed gradient is less harsh.

Expanding on the qualitative inspection, Figure 7 shows the reconstructions
for four images chosen to span a range of complexity (here we use the $\rm LL$
metric from Ulmo). Overall, the reconstructions are close to the original
images. Reflecting what we saw from Figure 6, as the masking ratio increases,
the model begins to miss or incorrectly guess some of the features. For
example, for image A and B for $t=75,p=75$ there are features which are
smoothed over or incorrectly reconstructed. In image C for $t=75,p=75$ in the
bottom right corner there is a small negative value (blue) feature that is
entirely missing in the reconstruction. When looking at the mask, we find that
the feature had been entirely masked out leaving Enki no information to
interpolate the feature. Despite this, Enki we are still able to recognize the
key patterns from the original image in the reconstruction indicating that
Enki is able to capture the overall structure of the cutouts. We continue to
observe these types of incorrect reconstructions in the $t=50$, $p=50$
examples, but to a lesser extent. In A we see for $t=50$, $p=5$ that the dot
of orange to the mid-left of the image is missing in the reconstruction, but
when looking at the blue gradients on the bottom-right and the swirl in the
middle, we find the features are not as smoothed out as they are in the
$t=75$, $p=75$ reconstructions.

Meanwhile, the images from D with LL=1110 are reconstructed well across all
models, and we suspect that complexity of a cutout will play a role in how
well Enki reconstructs images.

Figure 8:  The top image demonstrates how the Enki model trained at $t=75$
will frequently miss finer features in $p=75$ reconstructions. Notice how the
bottom right corner is almost entirely masked out. In the reconstruction, all
the structure of the original cutout is replaced with a smooth gradient. The
bottom image demonstrates how the $t=75$ model generalizes shapes, but looses
a lot of the finer details.  Figure 9: Gallery of different models
reconstructing cutouts at $p=10,30,50$. Notice that a model performs best when
reconstructing around the range for which it was trained. Models are still
able to capture the values of an image when reconstructing around 15-25% away
from where they were trained, but the patches become visible to the human eye.
This is especially noticeable in the $t=35$, $p=10,50$ reconstructions.

We next examined at how models reconstruct at masking ratios ($p$) outside of
their training ($t$). In Figure 9 we show an example of the $t=10,35,50,75$
models reconstructing at $p=10,30,50$. Here, we see issues arise with the
reconstructions. For the $t=10$ model, we see erros in the patches at $p=30$,
and at $p=50$ this effect is magnified. The $t=35$ model has this issue for
$p=10,50$, while the $p=30$ reconstruction does not. The $t=50$ model performs
significantly worse than the rest of the models at $p=10,30$, but does well at
$p=50$. This implies that there may be a certain offset between $t$ and $p$
for which a particular model of Enki will perform well.

Finally, the $t=75$ model does not have a single masking ratio within the
chosen range that produces acceptable results. At all levels of $p$ the
patches are visible, and there is an offset in values from the original image.

In their original paper, He et al. (2021) found that training at 75% masking
was ideal because it reduced redundancy when training, and forced the model to
learn more rather than extrapolating information from nearby patches. However,
for our purposes, we are more interested in accurate reconstructions over well
generalized ones. We have already seen that when a lot of smaller, finer
features are fully masked out, these features tend to either be smoothed over
or missing within the reconstructions. We demonstrate more in-depth examples
of this in Figure 8. In the top right corner of the upper image, Enki fails to
reconstruct the small band, and completely smooths over all the finer dynamics
in the bottom right. The bottom figure demonstrates another example of Enki
smoothing out finer details.

In addition, it was noted in the original MAE paper that different masking
techniques were tested. One of these techniques was a “block-wise” masking
strategy that masked out a large block of patches from the middle of the
image. While random patching is convenient for training due to the high
variability in masks, random patching does not accurately represent cloud
masks. This is especially true for the $p=10$ reconstructions where we have
many isolated patches scattered around the cutout. Block masking could more
accurately emulate large clusters of clouds blocking out large portions of the
image. We will consider this approach in future work.

Figure 10:  Average RMSE error for reconstruction of all masked $4\times
4\,\rm pixel^{2}$ patches at the location $i,j$ in the image. The $i$ and $j$
represent the position of the patch from the lower left corner. These are for
the validation set with $t=10$ and $p=20$ reconstructions. Clearly, patches at
the edge suffer significantly higher RMSE values. We attribute this to the
fact that edges and corners have less surrounding information compared to
patches in the center. For a real-world application of Enki, we would advise
ignoring the image-border for any reconstructions. All further calculations
and analysis presented in this thesis explicitly ignores the border.

###  4.2 Analysis of Reconstructions at the Individual Patch Level

Figure 11:  RMSE of patches vs the standard deviation of SSTa($\sigma_{T}$).
The black line is the one-to-one relationship where RMSE = $\sigma_{T}$. The
blue line is the RMSE of $t=10$, $p=20$ patches vs $\sigma_{T}$, and the red
line represents a (by-eye) fit to the blue line. We see that for
$\log_{10}\sigma_{T}>-2$ (0.02K), Enki’s RMSE is nearly a magnitude less than
the $\sigma_{T}$. This means the reconstructions are far more accurate than
random noise. At $\log_{10}\sigma_{T}<-2$, the RMSE flatten out. We speculate
that this is because spatial scale of the features become too small for
outside patches to assist in reconstruction.

An analysis of Enki’s performance was conducted at the individual patch level
to evaluate whether patch complexity or position affected Enki’s patch
reconstruction. We begin by calculating the RMSE of every patch of the
reconstructed $t=10$, $p=20$ images, and then taking the average of all
reconstructed patches at position $(i,j)$. As seen in Figure 10, the patches
along the border exhibit a systematically higher RMSE. This is especially true
at the corners which have the highest average RMSE. We speculate this is
because patches at the edges and especially the corners are reconstructing
with less information compared to patches towards the center. In essence, Enki
is forced to extrapolate instead of interpolate. A distinct example of this
can be see in Figure 6 where the gradient in the upper right corner is not as
pronounced as the one found in the original image. In some cases, when a
cutout has a small feature completely masked, the model will reconstruct the
image without the feature entirely.

While this issue may not be limited to edges and corners, Figure 6 indicates
that the edges and corners experience this issue more often than inner
patches. A blunt solution to handling these patches is to exclude them
altogether from the final cutout. While this throws out data, avoiding patches
more prone to this issue may be critical for real-world applications. For the
remainder of this thesis, we ignore the border in all further calculations and
analysis.

Continuing our analysis of the patches, we also examined Enki’s performance
based on the complexity of a patch. Specifically, we gauged the patch
complexity by measuring the standard deviation in temperature ($\sigma_{T}$).
We hypothesized that as patch complexity increases, we will also observe a
rise in reconstruction RMSE. After binning patches by $\sigma_{T}$, we
calculated the median reconstruction error (RMSE) of each bin (Figure 11).
Here, we plot the $\sigma_{T}$ of a patch against the RMSE of a patch on a
log-log scale.

The black dashed line represents the one-to-one relationship between
$\sigma_{T}$ and RMSE. If the results for Enki track this line it would mean
that the reconstructions are roughly equivalent to random noise. The blue line
is the RMSE of the $t=10$, $p=20$ patches plotted against the $\sigma_{T}$.
The red dotted line represents a by-eye fit.

We see that after a $\sigma_{T}$ of 0.02K ($\approx-2$ on the plot) that the
RMSE is smaller than the black line by nearly one magnitude. We do see that
there is a correlation between the RMSE and $\sigma_{T}$as they increase at a
similar rate, but the lower RMSE means that the constructions being produced
by Enki are excellent.

Before 0.02K, we observe that the plot flattens out at an RMSE of about
0.0012K (-3 on the plot). We suspect this is because as $\sigma_{T}$
decreases, the special scale of features in the SST patterns also decreases.
At the threshold of 0.02K, the spacial scale of the features become smaller
than the size of the patch, and the reconstructed field within the patch
becomes less reliant on information from outside the patch, resulting in patch
reconstructions that are pretty much random. This indicates that as we surpass
this threshold Enki’s reconstructions start using outside information when
reconstructing patches.

Figure 12: An example of the offset that can be found in the $t=75$ model at
$p=10$ and how adding a bias correction on the image can improve the final
reconstruction. We see that in the example without bias that the RMSE is
0.0267 and that there is an offset in this image towards 0.05K. In the 2nd
image, we subtract the bias of 0.0267 from the reconstructed patches, and see
an improvement in the RMSE down to 0.0162 and in the qualitative examination
of the reconstruction. Figure 13:  Histogram of the offsets in $t=75$
reconstructing $p=10$ images. The bias of the whole reconstruction set is
0.0267 which represented by the red line.  Figure 14: The bias that appears in
$t=10,35,50,75$ at $p=10,20,30,40,50$. We mark zero bias as a dotted red line.
In this plot, we see that in general, the $t=10$ model has a relatively low
bias across all Patch Mask Ratio. The next lowest is the $t=35$ which starts
out with a higher bias until it reaches 3 $p=30$ masking. The $t=75$ model
follows a similar pattern to the $t=35$ model. The $t=50$ stands out as having
a large bias at every $p$except $p=50$.

###  4.3 Bias Analysis

When examining models trained with a higher masking ratio ($t$) that are
reconstructing images of lower masking ratios (i.e. $t>p$), a systematic
offset is observed (e.g. Figure 13). This is most apparent when using the
$t=75$ model to reconstruct images with $p=10$ and for the $t=50$ model aside
from $p$=50. We calculate the offset in a single image as the median of the
absolute value of the difference between original and reconstructed pixels.
The mean offset of reconstructed set of images at a given $t$, $p$ is referred
to as the bias. For example, a single cutout reconstructed by the $t=75$ model
at $p=10$ may have an offset of 0.0131, but the bias of the whole
$t=75$,$p=10$ reconstruction set is 0.0269. When applying a correction to
reconstructions, we will use the bias because one cannot predict the correct
offset for a given cutout.

It is unclear why these offsets are present; we have yet to identify its
origin or even any strong correlation with image or patch properties. But when
the typical offsets are high enough, these lead to systematically large RMSE
values. We calculated the bias for all models at all percentage of masking
ratios, and plot them in Figure 14. In addition, we provide the values Table
2.

Figure 15: In this example, despite subtracting the bias (0.026 for
$t=75$,$p=10$) we see no improvement in the cutout. When we examine the RMSE,
we see that adding the bias correction made the cutout worse.

Figure 14 shows that the most prominent bias occurs for models with high
training mask ratios ($t>10$) applied to data with low masking ($p\leq 20$).
When examining $p=10$, a bias is present in all of the models aside from the
$t=10$ model. At $p=30$, we now see that the $t=35$ and $t=75$ model have a
much lower bias, nearing the performance of the $t=10$ model. We do observe
improvement in the $t=50$ model at $p=30$ as the bias increases by half from
$t=20$ to $t=30$.

As a whole, we note that the $t=10$ model’s bias always remains close to zero
while outside of the $p=50$.

Qualitatively and quantitatively, adding the bias improves the reconstructions
and lowers the RMSE. Occasionally, the offset of an individual image is much
larger than the bias, especially in images that have a larger $\sigma_{T}$. We
show an example of this in Figure 15 where despite adding a bias correction,
the RMSE is worse than if we had not added the bias. When calculating the
offset of this cutout, we find it to be 0.014. Despite this offset being half
of what the bias is, the change is barely visible (if at all).

Training Ratio | $p=10$ | $p=20$ | $p=30$ | $p=40$ | $p=50$
---|---|---|---|---|---
$t=10$ | 1.018e-0.6 | 1.214e-06 | 8.741e-0.6 | 8.027e-05 | 9.601e-05
$t=35$ | 0.0173 | 0.0077 | 0.0013 | -0.0013 | -0.0043
$t=50$ | -0.0288 | 0.0603 | 0.0298 | 0.00796 | -3.883e-6
$t=75$ | 0.0269 | 0.0102 | 0.00237 | -0.00183 | 0.00321
Table 2: The numerical values of all the biases. Notice how, for the
$t=10,35,50,75$ models, the greater $p$ is offset from the model’s training
ratio, the higher the bias. The only exception is the $t=10$ model, which has
a low bias for all masking ratios.

###  4.4 RMSE Analysis of Cutouts

We hypothesize that recontruction performance will degrade for images with
higher complexity, e.g., sharp fronts. To test this, we remove the edge
patches and apply a bias to all cutouts to ensure the best performance for all
models. For image complexity, we adopt the $\rm LL$ metric from Ulmo. For each
reconstruction set ($t$,$p$), we ‘binned the RMSE into 10 groups of ascending
LL and took the average of the RMSEs within those bins. We plot these results
in Figures 16 and 17.

When comparing the RMSE to the log likelihood ($\rm LL$) we see there is a
correlation between the two. As LL decreases, the RMSE increases. We also see
that for the the masking ratio has an effect on how well a model performs. For
$t=10$, as the masking ratio increases, the RMSE also increases.

The initial observation from Figure 17 is that the performance of the $t=50$
model was noticeably worse compared to the other models. This result matches
the bad reconstructions we have seen in Section 4.1 and the higher biases in
4.3. This disparity can be attributed to the lower learning rate used during
pre-training. We expect that increasing the number of training epochs or
relaunching pre-training with a higher learning rate will make the $t=50$
model more competitive with the other models. Despite its poor showings at
$p$, it is important to highlight that the $t=50$ model exhibited the best
performance when reconstructing images at $p=50$.

The $t=10$ model achieves the highest performance at $p=10$, followed by
$p=20$ and then $p=30$. In the case of the $t=35$ model, the best performance
is observed at $p=30$, followed by $p=40$ and $p=20$. The RMSE pattern of the
$t=75$ model follows a slightly different trend. It exhibits its worst
performance at $p=10$, shows improvement as it reaches $p=30$, and then
experiences a gradual increase in RMSE to $p=50$. Although the exact cause for
this behavior is uncertain, it is worth noting that the $t=75$ reconstructions
across all the examined masking ratios have had issues. This suggests that the
increasing masking ratio may contribute to the observed rise in RMSE after
$p=30$.

Once again matching what we observed in Figure 9, all models other than the
$t=10$ model exhibit their worst performance at $t=10$. These are also the
images that are observed to require the highest bias correction, and even with
such correction, their performance remains more than twice as poor as that of
the $t=10$ model. Clearly, these models have not learned the finer-scale
features apparent in the data.

We acknowledge that the $t=10$ model consistently has some of the lowest RMSE.
The only instances where we observe lower RMSE values than the $t=10$ model
are with the $t=35$ model at $p=40$ and the $t=50$ model at $p=50$. The $t=35$
model does, however, compare competitively with the $t=10$ model at $p=30$
masking with the RMSE appearing very similar.

While this would indicate that the $t=10$ model would be a good general model
to use for reconstructions, it is important to remember that when examining
the $t=10$ models qualitatively, we can see that the low RMSE values are
likely because it is able to estimate the values of the image well even at
higher masking ratios.

Overall Enki shows great promise in being able reconstruct SST images. While
some of the finer details can be lost especially in higher masking ratios, we
have seen countless examples in which Enki is still able to capture the
general structure of the cutout. In addition, for all masking ratios we
examined, there was always a model that had an RMSE lower than 0.05. We
believe that Enki shows great potential for reconstructing clouds masks in
real SST images.

Figure 16:  RMSE as a function of image complexity, here gauged by $\rm LL$ as
measured by the Ulmo algorithm with higher values indicating less complexity.
Irrespective of patch size, the average RMSE increase with image complexity
indicating the Enki model. Aside from the lowest bin in $\rm LL$ (most
complex), the RMSE of the reconstructed patches for $p\leq 30$ have average
RMSE $<0.04$ K which is comparable to the noise in remote sensing sensors like
VIIRS.  Figure 17:  Average RMSE for all of the image reconstruction sets. The
$t$=10 model performs well at all masking ratios $p$.  Figure 18:  Different
examples of how the $t=10$ model incorrectly reconstructs $p=40$ cutouts. In
the first, $t=10$ completely fails to guess the values of the original image
and fills in the patch with extreme values. In the second example, while
performing better than the first example, the patches are still visible. In
the third example, we see that $t=10$ has failed to both capture the values
and the structure of the patterns with patches have a blocky characteristic to
them.

##  5 Concluding Remarks

Our analysis reveals that Enki shows great promise as a method of
reconstructing cloud masked SST images. The qualitative examinations produce
sensible and usable results, and by exploring various training ratios, we aim
to identify the ratio that produces the most favorable and effective outcomes.

We find Enki is most accurate when handling reconstructions at lower masking
ratios with higher masking ratios running the risk of masking out too much
information which may lead Enki to smoothing over or missing masked features.
We also note, that patches at the edge of cutouts are more prone to error,
leading us to conclude that when reconstructing, it’s best to remove the edge
patches to reduce the chance of errors.

We do, however, acknowledge that there our some limitations in our current
analysis. The "clouds" we introduce into images are randomized, so they are
not a perfect representation of how real cloud masks will look. Going forward,
it would be good to test how Enki performs reconstructing LLC cutouts with
real cloud masks. Throughout this thesis, we also mainly focus on LLC data
which is free of a lot of issues that plague real data. Real remote sensing
data is noisy because of atmospheric noise and noise introduced by the
instruments, and it often runs the risk of unmasked clouds. The next step
would be to run Enki on real world data and see how it compares to older
methods of filling masked pixels like inpainting.

We also acknowledge that while Enki appears to be able to generalize well and
capture the main shapes of a dynamic, there is a lot of information that it is
unable to recreate accurately at a $p=75$ masking ratio. This becomes apparent
when comparing the RMSE of reconstructions at different masking percentages
where we see a steady increase in the RMSE as the masking ratio increases.
This means that while the models like the $t=75$ model can reconstruct up to
$p=75$ masking, we should be cautious because reconstructions at this level of
$p$ masking lose much of the details due to aggressive masking.

If we were to only use a single model, it would likely be the $t=10$ model.
Reconstructions would likely only go up to $p=20$, which if applied to VIIRS
data would increase the available cutouts by a magnitude (Figure 1). However,
since we know that different models perform best within a range of $p$, an
alternative approach is to incorporate multiple $t$ models within Enki and
select the appropriate model based on the cloud coverage of the cutout being
reconstructed. To determine the exact range where each model performs
optimally, further in-depth analysis is required.

Future work will encompass a range of tasks and objectives.

  * •

Our subsequent undertaking involves conducting analysis on VIIRS cutouts.
While the reconstructions of LLC cutouts displayed promising results, it is
crucial to determine if these outcomes can be effectively applied to real
data. This is because, despite the well-predicted SST patterns of VIIRS
determined by the LLC4320 model, it is important to acknowledge that the
LLC4320 data itself is generated and may have variations when compared to
real-world observations.

  * •

Moving forward, we also aim to conduct reconstructions utilizing real cloud
masks obtained from actual data. The current random masking, especially at
lower masking ratios, may not accurately represent the true cloud masking
scenario. Therefore, applying real cloud masks in our reconstructions has the
potential to yield different and more realistic results, providing valuable
insights into the performance of the reconstruction model.

  * •

Throughout this study, our training dataset and validation set consisted of
noise-free cutouts. While the primary objective of Enki is not to reconstruct
noise, we want to investigate whether the inclusion of noise has any impact on
the quality of reconstructions. By introducing noise into the input data, we
can assess the robustness and performance of Enki under more realistic
conditions.

  * •

Additionally, we consider conducting a follow-up study in which we set the
seed of the mask for reconstructions. This controlled experiment will enable
us to analyze the specific impact of each model on the reconstructed output,
without confounding factors introduced by random variations.

  * •

We may explore the use of different types of loss functions during training or
alternative equations to evaluate the performance of reconstructions. It is
worth considering that many images exhibit significant variations in stdT,
which means that the perception of reconstruction errors may differ based on
the magnitude of stdT. To address this, we could potentially evaluate a
model’s performance by calculating the RMSE normalized by the stdT of the
masked pixels in the original image. This approach would provide a more
contextually relevant assessment of the reconstruction quality, accounting for
the inherent variability in sea surface temperature across different images.

##  6 Acknowledgements

AA acknowledges J. Xavier Prochaska and P. C. Cornillon’s invaluable guidance
on this project, as well as David Reiman’s and Edwin Goh’s insightful
discussions. The authors acknowledge use of the Nautilus cloud computing
system which is supported by the following US National Science Foundation
(NSF) awards: CNS-1456638, CNS1730158, CNS-2100237, CNS-2120019, ACI-1540112,
ACI1541349, OAC-1826967, OAC-2112167.

## References

  * Baldi (2012) Baldi, P.: Autoencoders, unsupervised learning, and deep architectures, in: Proceedings of ICML workshop on unsupervised and transfer learning, pp. 37–49, JMLR Workshop and Conference Proceedings, 2012.
  * Bertalmio et al. (2001) Bertalmio, M., Bertozzi, A. L., and Sapiro, G.: Navier-stokes, fluid dynamics, and image and video inpainting, in: Proceedings of the 2001 IEEE Computer Society Conference on Computer Vision and Pattern Recognition. CVPR 2001, vol. 1, pp. I–I, 10.1109/CVPR.2001.990497, 2001.
  * Chelton et al. (2011) Chelton, D. B., Schlax, M. G., and Samelson, R. M.: Global observations of nonlinear mesoscale eddies, Progress in Oceanography, 91, 167–216, https://doi.org/10.1016/j.pocean.2011.01.002, 2011.
  * Dosovitskiy et al. (2020) Dosovitskiy, A., Beyer, L., Kolesnikov, A., Weissenborn, D., Zhai, X., Unterthiner, T., Dehghani, M., Minderer, M., Heigold, G., Gelly, S., Uszkoreit, J., and Houlsby, N.: An Image is Worth 16x16 Words: Transformers for Image Recognition at Scale, CoRR, abs/2010.11929, URL https://arxiv.org/abs/2010.11929, 2020.
  * Gallmeier et al. (2023) Gallmeier, K., Prochaska, J. X., Cornillon, P. C., Menemenlis, D., and Kelm, M.: An evaluation of the LLC4320 global ocean simulation based on the submesoscale structure of modeled sea surface temperature fields, arXiv e-prints, arXiv:2303.13949, 10.48550/arXiv.2303.13949, 2023.
  * He et al. (2021) He, K., Chen, X., Xie, S., Li, Y., Dollár, P., and Girshick, R. B.: Masked Autoencoders Are Scalable Vision Learners, CoRR, abs/2111.06377, URL https://arxiv.org/abs/2111.06377, 2021.
  * Isern-Fontanet et al. (2017) Isern-Fontanet, J., Ballabrera-Poy, J., Turiel, A., and García-Ladona, E.: Remote sensing of ocean surface currents: a review of what is being observed and what is being assimilated, Nonlinear Processes in Geophysics, 24, 613–643, 10.5194/npg-24-613-2017, 2017.
  * Jonasson and Ignatov (2019) Jonasson, O. and Ignatov, A.: Status of second VIIRS reanalysis (RAN2), in: Ocean Sensing and Monitoring XI, edited by Hou, W. W. and Arnone, R. A., vol. 11014, p. 110140O, International Society for Optics and Photonics, SPIE, 10.1117/12.2518908, 2019.
  * Levy et al. (2018) Levy, M., Franks, P., and Smith, K.: The role of submesoscale currents in structuring marine ecosystems, Nature Communications, 9, 10.1038/s41467-018-07059-3, 2018.
  * McWilliams (2017) McWilliams, J. C.: Submesoscale surface fronts and filaments: secondary circulation, buoyancy flux, and frontogenesis, Journal of Fluid Mechanics, 823, 391–432, 10.1017/jfm.2017.294, 2017.
  * Penven et al. (2006) Penven, P., Debreu, L., Marchesiello, P., and McWilliams, J. C.: Evaluation and application of the ROMS 1-way embedding procedure to the central california upwelling system, Ocean Modelling, 12, 157–187, https://doi.org/10.1016/j.ocemod.2005.05.002, 2006.
  * Prochaska et al. (2021) Prochaska, J. X., Cornillon, P. C., and Reiman, D. M.: Deep Learning of Sea Surface Temperature Patterns to Identify Ocean Extremes, Remote Sensing, 13, 744, 10.3390/rs13040744, 2021.
  * Prochaska et al. (2023) Prochaska, J. X., Beaulieu, C., and Giamalaki, K.: The rapid rise of severe marine heat wave systems, Environmental Research: Climate, 2, 021002, 10.1088/2752-5295/accd0e, 2023.
  * Prochaska et al. (2023) Prochaska, J. X., Guo, E., Cornillon, P. C., and Buckingham, C. E.: The Fundamental Patterns of Sea Surface Temperature, 2023.
  * Wu et al. (2017) Wu, F., Cornillon, P., Boussidi, B., and Guan, L.: Determining the Pixel-to-Pixel Uncertainty in Satellite-Derived SST Fields, Remote Sensing, 9, URL https://www.mdpi.com/2072-4292/9/9/877, 2017.