We investigate ocean circulation changes through the lens of data assimilation using a reduced-order model. Our primary interest lies in the Stommel box model which reveals itself to be one of the most practicable models that has the ability of reproducing the meridional overturning circulation. The Stommel box model has at most two regimes: TH (temperature driven circulation with sinking near the north pole) and SA (salinity driven with sinking near the equator). Currently, the meridional overturning is in the TH regime. Using box-averaged Met Office EN4 ocean temperature and salinity data, our goal is to provide a probability that a future regime change occurs and establish how this probability depends on the uncertainties in initial conditions, parameters and forcings. We will achieve this using data assimilation tools and DAPPER within the Stommel box model with fast oscillatory regimes.

MSC Classification: 37N10, 86A08, 65K99, 65P99

§ INTRODUCTION

Thermohaline circulation is the global current driven by density variations in the oceans. Its branch in the Atlantic ocean is known as the Atlantic Meridional Overturning Circulation (AMOC). In the northern part of the Atlantic it transports warm, salty water from the subtropics to the Arctic. Here it comes into contact with colder (sub)polar waters and sink forming the North-Atlantic Deep Water and flows southward in the deep ocean <cit.>. It plays a crucial role in the climate system, especially for the climate in northwestern Europe. Based on a simple two-box model, <cit.> suggested that AMOC can be described as a bistable system characterized by two distinct equilibrium states: TH (thermally-driven), representing the current circulation regime with sinking near the pole, and SA (salinity-driven), signifying a reversed circulation with equatorial sinking. The existence of multiple equilibria has since been confirmed by both reversal visible in the paleo-record as well as more sophisticated climate models, see e.g. <cit.>.

In the region where sinking occurs, temperature and salinity exert opposing effects on density. On one hand, heat exchange with the cold (sub)polar environment cools down the subtropical water that are being advected northward increasing its density. On the other, contact with melt water freshens it, decreasing its density. Currently, the former effect dominates. However, it is feared that climate change driven increases in temperature, ice melt and precipitation can tip the balance to the latter, thus reversing the circulation. The transition between the circulation regimes can occur through various mechanisms <cit.>. Beyond a critical point, parameter values may lead to bifurcation-induced tipping (B-tipping). Alternatively, the rate of parameter changes might trigger rate-induced tipping (R-tipping). Furthermore, system variability due to inherent noise could result in noise-induced tipping (N-tipping). Previous research, as documented in <cit.>, has established thresholds for tipping occurrence. In <cit.>, the Stommel box model is coupled to a sea-ice component in the polar box, and displays rate-induced tipping.

It is still an open question how close to such a tipping point the world is <cit.>. In 2023 alone, two studies on the topic <cit.> have received considerable attention in the press, but their results have been criticised <cit.>. The aim of this study is twofold. First, to explore whether it is possible to use observations of ocean temperature and salinity in combination with a reduced-order model in a process called data assimilation to estimate model parameters. Second, determine how the addition of forcing, both seasonal and climatological, to the Stommel model affects this probability of an AMOC reversal.

It is important to note that while one might assume that a two-box model would have limitations in fully describing tipping behavior, our approach, incorporating high-frequency oscillations, parameter variability, and data assimilation, brings us close to realistic scenarios. This methodology also sheds light on the feasibility of assimilating data to further understand the qualitative dynamics of our model. It offers a window of validity for our predictions, considering factors like the strength of the forcing. It's essential to clarify (and reiterate) that we do not assert that low-order models alone suffice for climate predictions. They serve as a means to grasp fundamental mechanisms and help capture essential physics.

The outline of the paper is as follows. First, we introduce the autonomous and non-autonomous Stommel model in Section <ref> and study the qualitative behavior of the solutions as time goes to infinity. In Section <ref>, we explain our data assimilation method and the assimilated observations. In Section <ref>, we conduct a twin experiment to investigate the potential of data assimilation to estimate model parameters. In Section <ref>, observations based on the Met Office Hadley Centre EN4 data set are assimilated in a non-autonomous Stommel model and the impact of the different forcings on the tipping probability is investigated. Finally in Section <ref>, we discuss our findings.

§ MATHEMATICAL ANALYSIS OF THE STOMMEL BOX MODEL: AUTONOMOUS AND NON-AUTONOMOUS PERSPECTIVES

§.§ The Stommel Box Model

The Stommel model comprises two fully mixed water containers. One represents the region of the North-Atlantic just north of the equator with temperature and salinity denoted as $(T_{e*}, S_{e*})$, while the other symbolizes the polar region with temperature and salinity as $(T_{p*}, S_{p*})$. (The asterisk subscripts here denote variables with a unit.) These containers are interconnected by both a capillary tube and an overflow area, facilitating water circulation while preserving a constant volume in each container. The relative salinities and temperatures of each box determine the flow of water, and when the flow of water reverses, we consider that to be a regime change.

A conceptual depiction of how Henry Stommel modeled thermohaline circulation. The circulation between the two connected bodies of well-mixed water in the Stommel model is driven by density gradients between the two boxes.

Initially, then, this scenario is described with four equations which are described below:

\begin{equation}
\begin{split}
\dx \dz \dy_{e} \frac{dT_{e*}}{dt_*} = & k_{T} \dx \dy_{e} (T_e^a - T_{e*}) + |\Psi_*|(T_{p*} - T_{e*})\\
\dx \dz \dy_{e} \frac{dS_{e*}}{dt_*} = & k_{S} \dx \dy_{e} (S_e^a - S_{e*}) + |\Psi_*|(S_{p*} - S_{e*}) \\
\dx \dz \dy_{p} \frac{dT_{p*}}{dt_*} = & k_{T} \dx \dy_{p} (T_p^a - T_{p*}) + |\Psi_*|(T_{e*} - T_{p*})\\
\dx \dz \dy_{p} \frac{dS_{p*}}{dt_*} = & k_{S} \dx \dy_{p} (S_p^a - S_{p*}) + |\Psi_*|(S_{e*} - S_{p*}) - m^{*} S_{p*}\\
\end{split}
\end{equation}

with $\dz$, $\dx$, the depth and the zonal width of the boxes and $\dy_{p}$, $\dy_{e}$ the meridional width of the  polar and equatorial box respectively, $k_{T}$ the temperature diffusion coefficient, and $k_{S}$ the salinity diffusion coefficient. The model incorporates the exchange of heat and salt within each box as a response to surface forcing, employing relaxation toward specified ocean surface temperature and salinity values $(T^a, S^a)$. Next to this, the model accomodates for the addition of fresh water supply to the polar box at a discharge rate of $m^{*}$.

The flow rate is given by the following equation:

\begin{equation}
\Psi_* = \gamma \dx \dz \frac{\rho_{p*} - \rho_{e*}}{\rho_0}
\label{eq:trans}
\end{equation}

where $\gamma$ is the advection coefficient and $\rho_0$ is a reference density. Since those are both constants, the flow rate is linearly related to the difference in water density between the two boxes. We note that if the flow is due to denser water in the polar box, as it is for the current TH regime, the flow rate is taken to be positive.

The Stommel system also assumes the following linear equation of state

\begin{equation} %equation of state
\rho_* = \rho_0(1 - \alpha_T(T_* - T_0) + \alpha_S(S_* - S_0))
\label{eq:eos}
\end{equation}

where a "0" subscript indicates reference value, $\alpha_T$ is the thermal expansion coefficient and $\alpha_S$ the haline coefficient.

But by instead looking only at the temperature and salinity relative to one another, we can reduce this down to a two variable system by subtracting the equations for the pole box from the equatorial box in the case that $m^{*}=0$:

\begin{equation}\label{eq:dimensionless}
\begin{split}
\frac{dT}{dt} =& \eta_1 - T(1+|T-S|)\\
\frac{dS}{dt} =& \eta_2 - S(\eta_3+|T-S|)\\
\Psi =& T - S
\end{split}
\end{equation}

with $T = \frac{\dz \gamma \alpha_{T}}{\overline{\dy} k_{T}} (T_{e*}-T_{p*})$,
$S = \frac{\dz \gamma \alpha_{S}}{\overline{\dy} k_{T}}(S_{e*}-S_{p*})$, $\Psi = \frac{1}{\overline{\dy} \dx k_{T}} \Psi_{*}$, $t=\frac{k_{T}}{\Delta z} t_{*}$ and $\overline{\dy}\overset{def}{=}\frac{\dy_{p} \dy_{e}}{\dy_{p}+\dy_{e}}$.

Here, $\eta_1,\eta_2$, and $\eta_3$ are parameters independent of the ocean state. Their full values are as follow:

\begin{align}
\eta_1 &= \frac{\dz \gamma \alpha_{T}}{\overline{\dy} k_{T}} (T^a_e-T^a_p) \label{eq:eta1}\\
\eta_2 &= \frac{k_S}{k_T} \frac{\dz \gamma \alpha_{S}}{\overline{\dy} k_{T}}(S^a_e-S^a_p) \label{eq:eta2}\\
\eta_3 &= \frac{k_S}{k_T}
\label{eq:eta3}
\end{align}

Setting $\frac{dT}{dt}=0$, $\frac{dS}{dt}=0$ in (<ref>), we get the following equilibrium solutions:

\begin{align*}
&& T = \frac{\eta_1}{1 + |\Psi|} && \text{and} && S = \frac{\eta_2}{\eta_3 + |\Psi|} \,, &&
\end{align*}

In the TH regime, we have $\Psi=|\Psi|$. In this case, the effect of the temperature difference dominates the effect of the salinity difference and the system is thermally driven. We find that

\begin{equation}\label{eq:Stom1}
\eta_2 = -\Psi^2 - \eta_3\Psi + \eta_1\left(\frac{\eta_3 + \Psi}{1 + \Psi}\right) \,.
\end{equation}

On the other hand, in the SA regime, we have $\Psi=-|\Psi|$. In this case, the transport due to the salinity difference dominates over the transport due to temperature difference. We find that

\begin{equation}\label{eq:Stom2}
\eta_2 = \Psi^2 - \eta_3\Psi + \eta_1\left(\frac{\eta_3 - \Psi}{1 - \Psi}\right)\,.
\end{equation}

The solutions from (<ref>) and (<ref>) for fixed $\eta_1$, $\eta_3$ can be shown together as a function of the dimensionless transport. From this, it is clear that for a range of $\eta_2$, three steady-state solutions for the transport are possible.

§.§ Lack of Periodic Solutions in the Stommel box model

In this section, we want to start by formally ruling out the existence of periodic solutions in the autonomous Stommel box model. Doing so will guarantee that in the non-autonomous model, when we do the twin experiments and begin assimilating real data, the only oscillatory behavior will come from the periodic forcing. The Poincaré-Bendixson theorem asserts that in a continuous dynamical system, any confined solution will inevitably approach either a fixed point or a limit cycle. In the next few paragraphs we will show that we get the former, not the latter in (<ref>), where we take $0<\eta_3<1$, $\eta_1>\frac{\eta_3}{1-\eta_3}$, with $T,S\geq0$.

To start, we have that

\begin{equation*}
\begin{split}
\frac{dT}{dt} =& f(T,S) = \eta_1 - T - T|T-S|\\
\frac{dS}{dt} =& g(T,S) = \eta_2 - \eta_3 S - S|T-S|
\end{split}
\end{equation*}

We define our Dulac function $\phi(T,S)=1$. Applying Dulac's criterion, we get:

\begin{equation}
\frac{\partial (\phi f)}{\partial T} + \frac{\partial (\phi g)}{\partial S} = \frac{\partial f}{\partial T} + \frac{\partial g}{\partial S}
= -1 - \eta_3 - 3|T-S|
\label{eq:dfdg}
\end{equation}

This is always strictly negative, and so by the Bendixson-Dulac theorem, we get that there cannot be any closed trajectories either entirely where $T>S>0$, as well as in the area where $S>T>0$. However, there could be a closed trajectory which crosses between the two areas, over the line $T=S$.

Now, assume we did have such a trajectory, a periodic solution $C$ which enclosed a region $D$ within $T,S\geq0$, and let us call the part part of $D$ where $T>S$, $D_1$, and the part where $S<T$ as $D_2$, so that $\partial D_1 = C_1 \cup C_3$ and $\partial D_2 = C_2 \cup (-C_3)$.

As <ref> allows for a continuous extension from $D_1$ ($D_2$) to its compact closure $\overline{D}_1$ ($\overline{D}_2$), $\int\int_{\overline{D}_1}(\frac{\partial f}{\partial T}+\frac{\partial g}{\partial S})\:dT\:dS$ ($\int\int_{\overline{D}_2}(\frac{\partial f}{\partial T}+\frac{\partial g}{\partial S})\:dT\:dS$) exists and is, as shown above, strictly smaller than zero. Applying Green's Theorem now gives:

\[\begin{split}
\int\int_{D}(\frac{\partial f}{\partial T}+\frac{\partial g}{\partial S})\:dT\:dS=&\int\int_{D_1}(\frac{\partial f}{\partial T}+\frac{\partial g}{\partial S})\:dT\:dS+\int\int_{D_2}(\frac{\partial f}{\partial T}+\frac{\partial g}{\partial S})\:dT\:dS\\
=&\int_{C_1\cup C_3}(-g\:dT+f\:dS) + \int_{C_2\cup (-C_3)}(-g\:dT+f\:dS)\\
=&\int_{C_1}(-g\:dT+f\:dS) + \int_{C_3}(-g\:dT+f\:dS)+\int_{C_2}(-g\:dT+f\:dS)\\
\end{split}\]

But rearranging this, the second and fourth terms simply cancel out. Furthermore, by definition, we have that $\frac{dT}{dt}=f(T,S)$ and $\frac{dS}{dt}=g(T,S)$, so we get:

\[\begin{split}
\int\int_{D}(\frac{\partial f}{\partial T}+\frac{\partial g}{\partial S})\:dT\:dS=&\int_{C_1}(-g\:dT+f\:dS) + \int_{C_2}(-g\:dT+f\:dS)\\
=&\int_C (-g\:dT+f\:dS)\\
=&\int_C (-g\frac{dT}{dt}+f\frac{dS}{dt})\:dt\\
=&\int_C (-\frac{dS}{dt}\frac{dT}{dt}+\frac{dT}{dt}\frac{dS}{dt})\:dt\\
=&\int_C 0\:dt\\
=& 0
\end{split}\]

But earlier we showed this was less than zero, so we have a contradiction, so no such loop $C$ can exist.

Thus, in the autonomous Stommel box model, it is impossible for a periodic solution to exist for $T,S \geq 0$.

§.§ The Non-Autonomous Stommel Model

We extend the Stommel model to include seasonal variations, represented by periodic changes in surface temperatures and salinities. As a consequence $\eta_{1}$, $\eta_{2}$ are no longer constant in time but oscillate. I.e.

\begin{align}
\eta_1 =& \eta_{1,auto} + B\sin(\Omega t) \label{eq:eta1_B}\\
\eta_2 =& \eta_{2,auto} + \hat{B}\sin(\Omega t) \label{eq:eta2_Bhat}
\end{align}

with $\eta_{\cdot,auto}$ represents the values used in the autonomous system in section <ref>.
Here, $B$ and $\hat{B}$ are distinct, but we define $A=B-\hat{B}$ for future use. Using that, our definition of $\Psi$, and the overall Stommel system, we can write:

\[\begin{split}
\frac{d}{dt}\Psi=&\frac{d}{dt}(T-S)\\
=&(\eta_{1,auto} + B\sin(\Omega t) - T(1+|T-S|))-(\eta_{2,auto} + \hat{B}\sin(\Omega t) - S(\eta_3+|T-S|))\\
=&\eta_{1,auto} - \eta_{2,auto} + \eta_3 S - T - (T-S)|T-S| + (B-\hat{B})\sin(\Omega t)\\
=&\eta_{1,auto} - \eta_{2,auto} + \eta_3 S - T - \Psi|\Psi| + A\sin(\Omega t)
\end{split}\]

then we substitute in $S=T-\Psi$, which ultimately gets us:

\begin{equation}\label{eq:nonaut}
\begin{split}
\frac{d}{dt}{\Psi}=&\eta_{1,auto} - \eta_{2,auto} + (T-\Psi)\eta_3 - T - \Psi|\Psi| + A\sin(\Omega t)\\
\frac{d}{dt}{T}=&\eta_{1,auto} + B\sin(\Omega t) - T(1+|\Psi|)
\end{split}
\end{equation}

for a system without climate change. In particular for a system in which $m^{*}=0$.

In order to get a full picture of the dynamics of (<ref>), it is thus fitting to analyze its $\Psi$-vs.-$\eta_{2,auto}$ bifurcation diagram. It's worth noting that the discontinuity surface $\Sigma_0 =  \{(T,S) \in \mathbb{R}^2_+ | T=S\}$ give rise to a Non-Smooth Fold (NSF) bifurcation, black star in <ref>, where the stable equilibrium branch terminates at $\eta_{2c}\equiv \eta_1\eta_3$. The NSF arises when the unstable equilibrium solution (dotted blue curve) and the focus (solid red curve) intersect with the discontinuity surface. This tipping point corresponds to the rapid transition from solutions near the temperature-dominated branch of nodes (solid blue curve) to the salinity-dominated branch of focus points and would arise particularly if we allow $\eta_2$ to be dynamic while having $A=B=0$. By dynamic, we mean slowly decreasing the bifurcation parameter $\eta_2$  of ”rate” $\varepsilon \ll 1$.

For the non-autonomous Stommel model ((<ref>)), for high frequency oscillatory forcing, there is an advance in the transition away from the branch of equilibria which depends linearly on the amplitude-to-frequency ratio $(A/\Omega)$ of the oscillations, thin yellow line in <ref>. The foundational work for this analysis can be found in a paper by Budd, Griffith, and Kuske <cit.>. Their study provides valuable insights into the expected behaviors of the parameters $\eta_1$ and $\eta_2$, as well as the equilibria, when incorporating seasonal variations. Building upon their findings, we extend the investigation by introducing data assimilation techniques to further enhance our understanding of these complex dynamics. By assimilating observational data into the model, we aim to explore how the combination of seasonal variations and data assimilation influences the system's behavior and equilibria.

A $\Psi$-vs.-$\eta_{2,auto}$ bifurcation diagram for the system in equation <ref> with $\eta_{1,auto}=3$, $\eta_3=0.1$. The solid red and blue lines are stable equilibria, while the dotted blue line represents an unstable equilibrium. They represent the equilibria for static parameters ($A=B=0$). The solid black circle indicates the system's saddle-node bifurcation, while the black star indicates the NSF bifurcation. The dotted black line emphasizes the discontinuity set $\Sigma_0$. The thin yellow line is an example solution (equilibria) to (<ref>) with $\Omega=100$, and $A=B=1$ (ratio $A/\Omega=0.01$).

§ INVERSION OF THE STOMMEL BOX MODEL

The solution of the Stommel box model in (<ref>) depends on the surface diffusion parameters $k_{T}$ and $k_{S}$ as well as the advection parameter $\gamma$. Given the simplified nature of the Stommel box model it is not possible to estimate appropriate values for these parameters from first principals. Data assimilation allows to reconstruct plausible values for these parameters from comparisons between model output and observations. In this section we will introduce the data assimilation methodology and the assimilated observations. The described model and methodology are implemented in the DA Python package DAPPER <cit.>.

§.§ Ensemble Kalman filter

Data assimilation (DA) is a procedure which combines a prior estimate of the state of a system, usually a forecast produced by a model, with actual measurements in order to produce a better prediction of the true state. Apart from the ocean state variables (temperature and salinity of the boxes in this case), the state also include the model parameters (diffusion and advection coefficients). This is referred to as the augmented approach. In this study we will use the ensemble transform Kalman filter (ETKF) <cit.> as data assimilation method. The ETKF is a variant of the Kalman filter (KF) <cit.>. In the KF it is assumed that the probability distribution at initial time is given by a Gaussian, that the model is linear and that it is combined with observations that possess a Gaussian error distribution. Under these assumptions the probability distributions remains Gaussian through time. Consequently, it can uniquely be described by a time-dependent mean and covariance. The KF is applied by cycling through the following two phases:

* Prediction Phase: Let $\Vec{x}' \in \mathbb{R}^{n}$ be the true state vector at the previous time step in a system with $n$ state variables. Since the true state is unknown, this vector is a random variable which we assume to have a Gaussian distribution with mean represented by the $n-$dimensional vector $\Vec{\mu}_x'$ and the $n\times n$ covariance matrix $\Sigma_{x'}$. If the model dynamics are linear, the probability distribution for $\Vec{x}$ at the current time time step is a Gaussian as well with mean and covariance given by $\Vec{\mu}_x = \mathbf{M} \Vec{\mu}_{x'}$ and $\Sigma_{x}=\mathbf{M} \Sigma_{x'} \mathbf{M}^{\rm{T}}+\mathbf{Q}$ with $\mathbf{M}$ the linear model and $\mathbf{Q}$ the covariance matrix for the additive model error (if any). The probability distribution of $\Vec{x}$ at the end of the prediction phase is referred to as the forecast.

* Update Phase: Let $\Vec{d}\in \mathbb{R}^{m}$ be a vector containing the $m$ real data observed at this time step. Due to errors in the measurements, this is also assumed to be a random variable with Gaussian distribution with covariance matrix $\Sigma_d$ and mean is $\Vec{\mu}_d=H\Vec{x}$, where $H$ is an observation matrix. It can be shown that the analysis distribution for $\Vec{x}$ given $\Vec{d}$ is a Gaussian too with mean and covariance given by:

\begin{align}
\Vec{\mu}_{\hat{x}} &= \Vec{\mu}_x + K(\Vec{d}-H\Vec{\mu}_x)\\
\Sigma_{\hat{x}} & = (I-K H)\Sigma_x = \Sigma_{x}
-\Sigma_{x} H^{\T} (H \Sigma_{x} H^{\T}+\Sigma_{d})^{-1} H \Sigma_{x}\, ,
\label{eq:kalman}
\end{align}

where $I$ is the identity matrix and $K$ is the so-called Kalman gain matrix. This matrix serves as a sort of weighting that determines how much the prediction $\Vec{x}$ is favored versus the observed data $\Vec{d}$. If the Kalman gain is close to zero, then $\Vec{d}$ has very little weight on the updated state vector $\Vec{\hat{x}}$. The Kalman gain can be calculated as such:

\[
K=\Sigma_x H^T(H \Sigma_x H^T + \Sigma_d)^{-1}
\]

The probability distribution for $\hat{x}$ is called the analysis.
By cycling between the two phases, one reduces overtime $\Sigma_{\hat{x}}$ and consequently can bring $\Vec{\mu}_{\hat{x}}$ reasonable close to the unknown true state $\Vec{x}$.

The Stommel box model as described in (<ref>) is however not linear and consequently straightforward application of the KF is not possible. To accomodate the non-linearity, the ETKF, a particular type of ensemble Kalman filter (EnKF) is used. EnKFs are Monte Carlo analogues of KF in which it the background probability distributions at each analysis time step is assumed to be Gaussian with mean and covariance estimated from an ensemble of model runs. This is enforced by defining the analysis ensemble as

\begin{equation}
\hat{X} = X+\Vec{1}^{\T} K(\Vec{d} - H\Vec{\mu}) + A Q \Lambda^{\frac{1}{2}} Q^{\T}
\end{equation}

with $\hat{X}$ a matrix having the analysis ensemble members as columns, $X$ the matrix with forecast ensemble members, $A=X(I-\Vec{1}^{\T}\Vec{\mu})$ the ensemble perturbations from the mean $\Vec{\mu}$ and $Q \Lambda Q^{\T}$ the singular value decomposition of $I - A^{\T} H^{\T} ((m-1)\Sigma_{d}+H A A^{\T} H^{\T})^{-1} H A$ <cit.>.

§.§ Log-normal transform

Temperature diffusion, salinity diffusion and advection coefficients are by definition positive. As the support of the Gaussian probability density function is unbounded, this is incompatible with the assumption that the probability density distribution of these parameters can be described by a Gaussian. In some of our trial runs with assimilating  realistic observations in the non-autonomous Stommel model, it was found that negative values were a common occurrence for these coefficients. To resolve the issue, the probability density distribution for these three parameters is assumed to follow a log-normal distribution. That is to say, the natural logarithms of the temperature diffusion, salinity diffusion and advection coefficients are assumed to come from a Gaussian distribution. Consequently, the ETKF described in <ref> can be used to correct these logarithms. Care should be exercised when obtaining a single best estimate of the coefficients from the ensemble. For a Gaussian distribution the most-likely estimate, or mode, of the distribution coincides with the mean. For the log-normal distribution this is not the case and the mode of the distribution of the coefficients is not given by $e^{\mu}$ but by $e^{\Vec{\mu}_{p}-\Sigma_{p} \Vec{1}}$ with $\Vec{\mu}_{p}$ the 3-dimensional vector containing the ensemble mean of the logarithms and $\Sigma_{p}$ the $3\times3$ matrix with their ensemble covariance <cit.>.

§.§ Observations and Model Parameters

Geometry of the two boxes in the model, box-average temperature and salinity, observation error variances and box-surface-averaged temperature and salinity necessary for the model's forcing are compiled based on the Met Office's EN4 quality controlled objective analysis <cit.> with bathymetric bias corrections <cit.>. The EN4 dataset is a dataset of several types of in-situ potential temperate and salinity measurements interpolated onto a 1x1 degree grid of the globe at 42 non-equally space depths on a monthly basis. From this dataset the vertical profiles in the North-Atlantic between $23.5 \, ^\circ \mathrm{N}$ and $80.0 \, ^\circ \mathrm{N}$ and between $80.0 \, ^\circ \mathrm{W}$ and $20.0 \, ^\circ \mathrm{E}$ that are at least 2 km deep are selected. These boundaries were chosen such model includes the North-Atlantic north of the intertropical convergence zone ($5-15 \,\unit{^\circ N}$). The precipitation in this zone leads to freshening of the ocean surface that the simple Stommel box model geometry cannot accommodate. Data from the period 2004-2022 are used. This selection is based on the fact that during this period extensive subsurface measurements were available thanks to the Argo program <cit.> as well as that they exclude a sharp spike of surface temperatures at the beginning of the century probably caused by processes our simple Stommel box model cannot represent.

With its two-box representation, the Stommel model is an oversimplification of the geometry and physics of the North-Atlantic basis. During experimentation it was found that results of the model and their physical plausibility depended strongly on the way the profiles were assigned to the different boxes. In order to ensure that the model is able to represent the physics of the thermohaline circulation, i.e. sinking of water in polar box driven by cooling, vertical profiles are divided between the boxes as follows. Those profiles in which surface temperature drops below the depth-averaged temperature for on average 1 month per year or more are assigned to the polar box, the other profiles to the equator box.

Polar     Equator

$\dx$ $\mathrm{[km]}$     5210     5210

$\dy$ $\mathrm{[km]}$     166     3853

$\dz$ $\mathrm{[m]}$     4280     4280

The zonal width, meridional width and depth of the two boxes used in the Stommel model based the EN4 data.

Volume $V$ of each of the boxes is calculated by summing the volume of all grid cells in the box. Similarly, surface area $A$ is calculated by summing the area of the EN4 grid cells in the top layer in both boxes. The depth of both boxes is then calculated as $\Delta z = \frac{V_{p}+V_{e}}{A_{p}+A_{e}}$. The cross-section $A_{c}$ between boxes is derived by summing the vertical area between profiles assigned to the polar box and those assigned to the equatorial box. An example assignment is shown in Figure <ref> From this zonal and meridional size of the boxes is defined as $\Delta x = \frac{A_{c}}{\Delta z}$ and $\Delta y_{p}=\frac{A_{p}}{\Delta x}$, $\Delta y_{e}=\frac{A_{e}}{\Delta x}$. The values resulting from this exercise are shown in <ref>.

An example of EN4 profile assignment in the year 2010. Orange points signify profiles assigned to the equatorial box, and blue points signify profiles assigned to the polar box.

Time-averaged obs.     Obs. error std. dev.

$T_{p*}$ $[\celsius]$     3.0      0.3

$T_{e*}$ $[\celsius]$     9.4     0.4

$S_{p*}$ $[\ppt]$      34.86     0.05

$S_{e*}$ $[\ppt]$     35.14     0.06

Time-average of the assimilated temperature and salinity observations together with the observational standard deviation, i.e. the diagonal of $\Sigma_{d}^{\frac{1}{2}}$.

For each month $t$ and each of the two boxes the box-averaged temperature is calculated as

\begin{equation}
T_{\cdot *}(t) = \big( \sum_{i} T_{i*}(t) \frac{V_{i}}{\sigma_{T,i*}^2(t)} \big)
\big( \sum_{i} \frac{V_{i}}{\sigma_{T,i*}^2(t)} \big)^{-1}
\label{eq:T_mean}
\end{equation}

where $i$ runs over all EN4 grid cells in the box below the surface layer, $V_{i}$ is the volume of the grid cell, $\sigma_{T,i}$ is the uncertainty in temperature as specified in the EN4 data set. The variance in the value is calculated as

\begin{equation}
\sigma_{T*}^{2}(t) = \big( \sum_{i} V_{i} \big)
\big( \sum_{i} \frac{V_{i}}{\sigma_{T,i*}^2(t)} \big)^{-1}
\label{eq:T_sigma}
\end{equation}

These calculations are repeated for salinity $S$.  The values after January 2004 are used to create the observations to be assimilated in $\Vec{d}(t)=[T_{p*}(t),T_{e*}(t),S_{p*}(t),S_{e*}(t)]^{\mathrm{T}}$. The $4 \times 4$ observational error covariance matrix $\Sigma_{d}$ is assumed to be diagonal. The values on the diagonal are chosen to be the maximal values for $\sigma_{T*}^{2}(t)$ and $\sigma_{S*}^{2}(t)$ in each box. The resulting time-average values for the observations as well as the observational error standard deviation found this way are shown in <ref>.

Field     $\cdot_{0}$     $\cdot_{\cos}$     $\cdot_{\sin}$     seasonal amplitude

$T^{a}_{p}$ $[\celsius]$     2.1     -2.0     -1.3     2.4

$T^{a}_{e}$ $[\celsius]$     16.2     -2.0     -2.0     2.8

$S^{a}_{p}$ $[\ppt]$     33.84     0.33     0.27     0.43

$S^{a}_{e}$ $[\ppt]$     36.90     -0.01     -0.03     0.03

Regression coefficients found by fitting model <ref> to the surface observations for temperature $^{a}$ and salinity $^{a}$ in polar $\cdot_{p}$ and $\cdot_{e}$ boxes. Also shown is the seasonal amplitude $\sqrt{\cdot_{\cos}^2+\cdot_{\sin}^2}$

The calculations in <ref> and <ref> are repeated but now with $i$ running over all grid cells in the surface layer of a box. The signal is then represented as

\begin{equation}
T^{a}(t) \approx T^{a}_{0} + T^{a}_{\sin} \sin(\frac{2\pi t}{\tau}) + T^{a}_{\cos} \cos(\frac{2\pi t}{\tau})
\label{eq:T_regression}
\end{equation}

with $\tau=1 \, \mathrm{year}$ and $T_{0}$, $T_{\cos}$, $T_{\sin}$ determined by a weighted linear regression against the box-averaged surface values using $\sigma_{T*}^{-2}(t)$ as weights. The regression coefficient for surface temperature and their equivalents for surface salinity are given in <ref>.

$\Omega$     $B$     $\hat{B}$     $A$

dimensional     $2\pi$ $\mathrm{rad\,yr^{-1}}$     0.66 $\celsius$     0.46 $\ppt$     -

scale     $\frac{k_{T}}{\dz}$     $\frac{k_{T}}{\alpha_{T} \gamma \dz}\frac{\dy_{p} \dy_{e}}{\dy_{p} + \Delta y_{e}}$
$\frac{k_{T}}{k_{S}}\frac{k_{T}}{\alpha_{S} \gamma \dz}\frac{\dy_{p} \dy_{e}}{\dy_{p} + \dy_{e}}$     -

nondimensional     $950$     $1.98 \cdot 10^{-1}$     $1.51 \cdot 10^{-1}$      $4.6 \cdot 10^{-2}$

Relation between the coefficients $A$, $B$, $\hat{B}$ in <cit.> and the surface forcing obtained from the EN4 observation data.

The regression coefficients in <ref> can be related to the perturbations added to $\eta_{1}$ and $\eta_{2}$ in Equation (<ref>) and Equation (<ref>). These relations can be found in <ref>.

Parameters, forcings, initial conditions and their sources.

Data     Values     Source     Details
Temperature (ocean temperature)        Pole: $T^{*}_{p}=3.15 \,\unit{\celsius}$, Equator: $T^{*}_{e}=9.42 \,\unit{\celsius}$
Initial values from EN4 data

Salt (ocean salinity)        Pole: $S^{*}_{p}=34.84 \,\unit{ppt}$
Equator: $S^{*}_{e}=35.21 \,\unit{ppt}$
Initial values from EN4 data

temp_diff (surface heat mixing coefficient)        $k_{T}=\frac{\kappa_T}{z_{mixing}}$ with $\kappa_T=10^{-4}\mathrm{m^2 s^{-1}}$

salt_diff (surface salinity mixing coefficient)        $k_{S}=\eta_{3} k_{T}$
<cit.>     $\eta_3=0.1$

$z_{mixing}$ (mixing layer depth)        $10 \,\unit{m}$
Based on 2019 Hadley data

$\dx$ (width of polar and equatorial boxes in lateral direction)         $5210 \,\unit{km}$
Calculated from EN4 dataset; see Table <ref>

$\dy$ (width of polar and equatorial boxes in longitudinal direction)     Pole: $\dy_{p}=166 \,\unit{km}$, Equator: $\dy_{e}=3853 \,\unit{km}$
Calculated from EN4 dataset; see Table <ref>

$\dz$ (depth of polar and equatorial boxes basin)       $4280 \,\unit{m}$
Calculated from EN4 dataset; see Table <ref>

$\rho_0$ (reference density)       $1027 \,\unit{kgm^{-3}}$

$S_0$ (reference salinity)       $35 \,\unit{ppt}$

$T_0$ (reference temperature)     $10 \,\unit{\celsius}$

$\alpha_T$ (thermal expansion coefficient)     $\frac{0.15 \,\unit{\celsius^{-1}}}{\rho_{

$\alpha_S$ (haline expansion coefficient)     $\frac{0.78 \,\unit{ppt^{-1}}}{\rho_{

Initial estimate $Q_{overturning}$ (meridional overturning flux))     $18 \,\unit{Sv}$

$\gamma$ (advective transport coefficient)     $6.9 \cdot 10^{-2} \,\unit{ms^{-1}}$
Chosen such that $Q_{overturning}=18 \,\unit{Sv}$ at $t=0$.

$T^a(t)$     Oscillates over time
Calculated from EN4 dataset;  see table <ref>

$S^a(t)$     Oscillates over time
Calculated from EN4 dataset; see table <ref>

V_ice     $2.9 \cdot 10^6 \,\unit{km^3}$

Greenland ice sheet melt period $t_{melt}$     $10000 \,\unit{y}$     <cit.>

Equatorial warming rate     $0.03 \,\unit{\celsius y^{-1}}$      <cit.>

Polar warming rate
$0.06 \,\unit{\celsius y^{-1}}$      <cit.>

Model time step     1 month

Initial conditions for the ocean temperatures and salinities in the different ensemble members are drawn from a Gaussian distribution having the EN4 box-averaged temperatures and salinities for January 2004 as mean and $\Sigma_d$ as covariance. Initial values for the logarithm of the diffusion and advection parameters are chosen such that the mode of the distribution corresponds to $k_{T}=1.0 \cdot 10^{-5}  \unit{\frac{m}{s}}$, $k_{S}=1.0 \cdot 10^{-6} \unit{\frac{m}{s}}$ and $\gamma=1.3 \unit{\frac{m}{s}}$. The logarithms are perturbed by adding realizations from a Gaussian with a standard deviation of $0.26$. This corresponds to the assumption that the relative error in these model parameters is $30\%$.  Values for the other parameters used in the model can be found in Table <ref>. The numerical scheme for Stommel model uses a finite volume approach with an explicit 4th-order Runge-Kutta scheme with a time step of 1 month.

§ TWIN EXPERIMENTS

In this section we will evaluate whether the data assimilation system is capable of correcting ocean temperatures and salinities as well as reconstructing the Stommel model parameters. In addition to this we want to see what happens with the AMOC circulation if the effects of climate change are added to the model. This is done by running an Observing System Simulation Experiment (OSSE), or twin experiment, in which one of the ensemble members is set apart as“truth” and withheld from the DA system. Artificial observations are created from the “truth” by adding perturbations drawn from a Gaussian distribution with covariance $\Sigma_{d}$ with these perturbations mimicking observation errors. After running the experiment, the forecasts and analyses can be compared with the “truth” to quantify performance of the DA system.

No climate change, no DA, temperature difference

No climate change, no DA, salinity difference

No climate change, no DA, phase portrait

No climate change, synthetic data, temperature

No warming, synthetic data, salinity

No climate change, synthetic data, phase portrait

Climate change, no DA, temperature difference

Climate change, no DA, salinity difference

Warming, no DA, phase portrait

Climate change, synthetic data, temperature

Climate change, synthetic data, salinity

Climate change, synthetic data, phase portrait

Forecast time series showing temperature and salinity differences between equatorial and pole boxes for 100 ensemble members. Black paths: synthetic truth, orange paths: ensemble mode. Ensemble trajectories plotted in dimensionless temperature vs salinity phase portraits (initial conditions: green dots). Orange dots with yellow error bars: assimilated observations with error standard deviation.

We conduct simulations under four distinct scenarios. In the first scenario initial conditions and model parameters are perturbed in the different ensemble members but no data assimilation is used. The second scenario mirrors the first but now ocean temperatures and salinities are assimilated each month during the period 2004-2022. For the third scenario the effects of global warming are taken into account. In this study we take two processes related to climate change into account. First, that the average temperature increases faster near the pole  ($0.06 \unit{\frac{\celsius}{year}}$ for polar box) than at lower latitudes ($0.03 \unit{\frac{\celsius}{year}}$ for equatorial box). Second, that the Greenland ice sheets melts add a constant rate $m^{*}=\frac{V_{ice}}{t_{ice}}$ with $V_{ice}$ and $t_{ice}$ given in Table <ref>.  In order not to interfere with the DA, warming and melting only start after 2022, i.e. after the end of the DA period, and is applied uniformly to all ensemble members. The fourth scenario matches the third and includes data assimilation. The model time step is 1 month and DA is carried out monthly. <ref> shows the model ocean temperatures and salinities over 100 years for 100 ensemble members together with its most-likely value (the ensemble mean) and the artificial truth. Shown ensemble members are forecast solutions, i.e. solution just before DA is applied if it is used. Temperatures and salinities in the 3rd column are non-dimensionalized using the scaling constants in Section <ref>. Blue curves indicate the thermohaline (TH) regime, while red curves represent the salinity-driven (SA) regime. As the direction of the AMOC circulation and the equilibrium solutions of the Stommel model depend solely on the differences between the boxes, only the differences between the equatorial and pole box are shown.

As seen in Figures <ref>-<ref>, none of the ensemble members predict a change in flow direction in the experiment with no climate change and no data assimilation. This is notwithstanding the large uncertainties in the model parameters. One can see that there is a significant spread between trajectories in the time plots, but the ensemble remains monomodal, i.e. the ensemble members do not form different clusters. This is a necessary condition for the EnKF to function properly. Figures <ref>-<ref> show that assimilating temperatures and salinities for the period 2004-2022 reduces the spread bringing all ensemble members closer to the truth. Trajectories become and stay tightly connected in the time and phase plots. This indicates that DA is successful in correcting the temperatures and salinities.
When the effects of climate change are applied in Figures <ref>-<ref> the solutions no longer converges towards an equilibrium value. Instead temperature and salinity keep increasing reflecting the introduced changes in forcing. Apart from this, the results are qualitatively the same as in the first two scenario's: the circulation remains in TH mode and the DA reduces the errors in the ensemble members.

No warming temp. flux

No warming salinity flux

No warming advective flux

The temperature, salinity, and advective flux coefficients of the ensemble plotted over the assimilation period for an experiment assimilating 18 years of synthetic data without any global warming effects (the experiment with global warming is effectively identical). The red curve represents the ensemble mode, while the green line represents the true coefficient value.

Polar Temp. RMSE

Polar Temp. Std. Dev.

Polar Salinity RMSE

Polar Salinity Std. Dev.

Equatorial Temp. RMSE

Equatorial Temp. Std. Dev.

Equatorial Salinity RMSE

Equatorial Salinity Std. Dev.

Temp. Diffusion Coeff. RMSE

Temp. Diffusion Coeff. Std. Dev.

Salt Diffusion Coeff. RMSE

Salt Diffusion Coeff. Std. Dev.

Advective Coeff. RMSE

Advective Coeff. Std. Dev.

RMSE and std. dev. of assimilated quantities in twin experiments. Black dotted line: end of assimilation. Colors: blue - no climate change, no assimilation; orange - climate change, no assimilation; green - no climate change, assimilation; red - climate change, assimilation. Green and blue may overlap with red and orange, appearing absent in some plots.

To quantify the behavior of the ensembles in <ref> in more detail, the ensemble standard deviation, or ensemble spread, and RMSE for the different experiments are shown in <ref> for the different state variables. The first thing that stands out is that ensemble spread and RMSE grows towards the end of the simulation period, apart from those for the equatorial temperature (<ref>f).  Here, even in the absence of DA, the spread decreases. This in contrast to the temperature in the pole box which nearly doubles (<ref>b). The explanation for this can be found in the large ratio of the box areas (see Table <ref>): the surface area of the equatorial box is 23 times as large as that of the pole box. Consequently, the surface heat flux dominates the advective heat flux in this box resulting in a rapid equilibration towards the surface temperature which is the same for all ensemble members. In the smaller pole box equilibration to the surface temperature is counteracted by the uncertainty the advection parameter and the uncertainty in the advective flux resulting from that. The second noticeable point is that DA is successfull in reducing the error in the model parameters. Errors for the advection parameter (<ref>m) and temperature diffusion parameter (<ref>i) are decimated while that of the salinity parameter is approximately halved (<ref>k). This reflects that salinity diffusion parameters is smaller. Consequently, it takes the system more time to equilibriate with surface salinity forcing and consequently the DA is less sensitive to changes in the salinity diffusion coefficient than in the advection and temperature diffusion coefficient. The fact that by the end of the DA period the spread in the salinity diffusion parameter is still decreasing (<ref>l), whilst that in the temperature diffusion parameter (<ref>j) and advection parameter ((<ref>n) has leveled off suggest that further RMSE reduction for the former could be achieved by extending the DA period. Finally, for a properly calibrated ensemble Kalman filter the expectation value of the ensemble spread and root mean-square error (RMSE) should be equal. Comparison of the RMSEs and ensemble standard deviations at the end of the DA period shows they are of similar magnitude. Thus confirming that DA is performing as expected.

Temperature and salinity time series, as well as a phase portrait for an experiment with yearly warming of $0.07 \,\unit{\celsius}$ per year in the equator ($0.14 \,\unit{\celsius}$ per year at the pole) and a melt period of 1000 years without DA. Trajectories turning red indicate a predicted flipping of flow direction. The black trajectories represent the synthetic truth values, and the orange represents the ensemble mode.

All ensemble members in <ref> are in the TH circulation mode. In order to see whether the setup is capable of producing solutions with the SA circulation, the impacts of the climate change are scaled up. More precisely, the warming rate in the equatorial box is increased to $0.07 \,\unit{\celsius}$ per year (the polar box warming rate is raised to $0.14 \,\unit{\celsius}$ per year) and the melting period for the Greenland ice sheet is reduced from 10,000 to 1,000 years. Once again, these climate change driven effects are only switched on at the end of the DA period otherwise the DA system would try to correct errors in the forcing by modifying the model parameters. In <ref> temperature (salinity) between-box difference in all ensemble members initially drops and then decreases (increases) due to the effects of climate change. However, <ref> shows that depending on initial conditions and the model parameters of the ensemble member, there is a rapid increase in the temperature and salinity between-box difference eventually resulting in a reversal of the AMOC circulation direction. So, the setup is capable of producing solutions with a circulation change.

To summarize, in this section we tested the ability of our setup to predict the truth using a twin experiment. We find that the setup produces an ensemble that covers the truth and can, with the sufficient forcing, produce solutions both with the TH as well as the SA circulation. DA is successful in reducing the spread and errors in both temperatures, salinities, diffusivity parameters and advection parameter. The 18 year DA period used is sufficient to reduce the uncertainty in temperature diffusivity and advection parameter to the lowest value attainable with this setup, but the uncertainty in the salinity diffusivity could have been reduced further by extending the DA period. Having established that our DA setup has the ability to recover the model parameters from temperature and salinity observations, we will use those observations to find realistic values for these parameters in Section <ref>.

§ ASSIMILATING BOX-AVERAGED OBSERVATIONS IN THE NON-AUTONOMOUS STOMMEL MODEL

In this section, we run simulations using (<ref>) while assimilating box-average observations compiled from the EN4 dataset as described in Section <ref>. In this section we pursue a two-fold aim. First, we attempt to find realistic values for the salinity diffusivity $k_S$, temperature diffusivity $k_T$ and advection parameter $\gamma$ as well as the non-dimensional values $\eta_{1}$, $\eta_{2}$ and $\eta_{3}$ which depend on these parameters. Second, having found realistic values, we continue running the box model forward in time to determine the likelyhood of a circulation reversal happening in this century both for the case with climate change and the case without.

§.§ Interfered Model Parameters

The (a) temperature diffusivity, (b) salinity diffusivity and (c) the advection parameter for the different ensemble members as function of time.  The red curves represents the ensemble mode.

Parameter      Most-likely     90%-confidence interval

$k_{T}$  $[10^{-7} \mathrm{ms^{-1}}]$      $9.0$      $[8.5, 9.5]$

$k_{S}$  $[10^{-7} \mathrm{ms^{-1}}]$      $1.9$      $[1.6, 2.2]$

$\gamma$ $[10^{-2} \mathrm{ms^{-1}}]$     $6.9$     $[6.6, 7.1]$

Most-likely estimates of model parameters before and after assimilating monthly boxed-averaged temperatures and salinities for the period 2004-2022 together with their $90$%-confidence interval.

The three different model parameters are shown in <ref> as function of time. As the effects of climate change are never applied during the DA period (see <ref>), the figure will look the same for both experiments with and without climate change. Notice that we are still assuming that within the conceptual constraints of the Stommel box model the true value of these parameters is fixed in time. The time dependence reflects that the probability distribution for the truth is changing as more and more observations get assimilated. Had we used the Rauch-Tung-Striebel smoother <cit.> we would have been able to correct the parameters in the beginning of the DA window using observations at the end of the window. However, as this would have not changed the final estimate, i.e. the value for the model parameters at the end of the DA window used to propagate the ensemble members forward, we have abstained from the additional effort. The figure shows an initial overestimation of the value of all three coefficients: the mode of the ensemble and all ensemble members exhibit a decreasing trend that starts to level off near the end of the DA period. Similar to the twin experiment in <ref> the salinity diffusion parameter has not fully converged near the end of the DA period. The ensembles for all three parameters show a steadily decreasing spread indicating that the uncertainty of the parameter's true values decreases. The final estimates for the parameters, i.e. their values at the end of the DA period, together with their 90% confidence interval based on the ensemble spread are given in <ref>.

The values of parameters $\eta_{1}$ (Equation (<ref>)), $\eta_{2}$ (Equation (<ref>)) $\eta_{3}$ (Equation (<ref>)) as function of time for the experiment with climate change.

The values for the non-dimensional parameters $\eta_{1}$, $\eta_{2}$ and $\eta_{3}$ as appearing in Equation <ref>
are shown in <ref>. The annual oscillations clearly visible in $\eta_{1}$ and $\eta_{2}$ are the result of seasonal cycle imposed on the surface forcing as Equation (<ref>). As climate change warms the polar box faster than the equatorial one, the temperature difference between the two decreases over time resulting in a decreasing slope in $\eta_{1}$ after the DA period. The annually averaged values for the ensemble mode of $\eta_{1}$, $\eta_{2}$ and $\eta_{3}$ at the end of the DA period are $4$, $0.8$ and $0.2$ respectively. These are in agreement with typical values explored in other Stommel box model literature <cit.>.

§ FUTURE CIRCULATION

Forecast time series starting in 2004 representing the equator-pole difference in (a) temperature and (b) salinity for the different ensemble members as well as (c) the ensemble trajectories in a dimensionless temperature vs salinity phase portrait. The orange dots with bright yellow errorbars represent the observed values and their standard deviation. Solution with a TH circulation are shown in blue, those with a SA circulation (if any) in red.

The differences between the forecast temperature (forecast salinity) in the equatorial box and the polar box for the 100 ensemble members for the experiment without climate change are shown in <ref>. It can be seen that during the DA period (first 18 years) both temperature and salinity differences follow the trend in the observations. However, they fail to reproduce the multiyear oscillation that is visible in the salinity and especially in the temperature observations. This is because such multiyear oscillations are missing in the surface forcing. This results in the model overestimating the salinity difference between equatorial and polar box thus underestimating the transport between the boxes. As explained in <ref>, the autonomous solutions converge to one of at most two possible stable equilibria each associated with a different circulation regime. In the non-autonomous system, these equilibrium values will oscillate in time and the solution always “chases” one of them or switches from chasing one to chasing the other. The figure shows that in the absence of climate change the circulation remains in the TH regime for all ensemble members until at least 2104. Furthermore, the solutions tend to values that are close to current conditions: in the Stommel box model the equilibrium solutions for the temperature differences lie below the difference at the end of the DA period. However, the differences are still within the range of contemporary temperature differences. The equilibrium solutions for salinity differences are larger than differences at end of the DA period, partially offsetting the decreasing trend over the DA period. The relative spread at the end of the century is also larger for the salinity differences than for the temperature differences. This can be explained by our observation in <ref> that the spread in salinity diffusion parameter is larger than that in the temperature diffusivity parameter. The equilibrium values depend on the model parameters and consequently the uncertainty for future ocean salinities is relatively larger than that in ocean temperatures.

As <ref> but now with effects of climate change (see text).

The solutions in which after the DA period the polar surface temperature increases with $0.06 \,\unit{\celsius y^{-1}}$, the equatorial surface temperature with $0.03 \,\unit{\celsius y^{-1}}$ and in which the Green land ice sheet melts with a constant rate over 10 thousand years are shown in <ref>. Compared to <ref> the salinity differences increase after the DA period (<ref>b). This is due to the ice melt. Next to this, the decrease in temperature difference is smaller due to the differential increase in the surface temperature between the two boxes (<ref>a). Apart from these points, there are no qualitative differences between the figures. In particular, none of the ensemble members changes circulation regime.

Percentage of ensemble members that undergo circulation regime change. All runs use a similar setup to the experiment with climate change and real data assimilation, but using different ice melt periods (x-axis) and different warmings rates for the equator box (y-axis). The warming rate in the polar box is twice that of equator box.   Bright red denotes mostly flipping, while blue represents mostly remaining in the starting TH regime.

The actual future rate of surface warming and ice melt will depend on, among other things, future carbon emission and are consequently uncertain. Therefore the DA experiment has been repeated several time for a range of plausible ice melting and between-box surface warming rate differences <cit.>. The percentage of ensemble members in which the circulation regime changes from TH to SA are shown in <ref>. The figure shows that after DA the uncertainty in the regime is small: either all ensemble members are in the TH regime or in the SA regime. Next to this, the close-to-vertical slope of the boundary between the two regimes indicates that faster ice melting is more influential in pushing the solutions to the SA regime than faster warming over the pole than equator. Measured melt rates over the last decade vary between 286 to 487 $\unit{Gt y^{-1}}$ which corresponds to a melt period of 7-12 thousand years <cit.>. This would place the future climate firmly within the zone of <ref> in which circulation the AMOC does not reverse direction.

The value of the between-box transport as function of time for the experiment with climate change. The vertical dashed line marks the end of the assimilation period. The horizontal dashed line indicates observed current transport <cit.>.

The transport between the two boxes is shown in <ref>. Transports are an order of magnitude smaller than the typical observational transect estimates over this period that lie in the range $10--25 \,\unit{Sv}$ <cit.>. The decrease is however not the result of climate change, but occurs during the DA period. I.e. the DA is not capable of simultaneously correcting temperatures, salinities and transports. The source of this discrepancy might lie in our box definition. The southern boundary of the equatorial box is placed at $23.5 \unit{^\circ S}$. Moving the border farther south would increase between-box temperature difference and thus produce higher transports even when using the same model parameters. The underestimation of current transports implies that the model is closer to the tipping point than the real AMOC. Consequently, the ice melt and warming rates in <ref> at which the circulation flips are in all likelihood an overestimation.

§ CONCLUSION AND DISCUSSIONS

We have applied the Stommel 2-box model to replicate the Atlantic meridional overturning circulation and investigated its solutions, the feasibility of retrieving the temperature diffusivity, salinity diffusivity and advection parameter by assimilating ocean temperature and salinity parameters and to forecast the likelihood of circulation regime change in the 21st century.
First, we showed using the Poincaré-Bendixson theorem that the autonomous model does not have any periodic solutions. Second, we introduced seasonal variability and two of the effects of climate change into the model: warming of the ocean surface at different rate for the pole and the equator and the melt of the Greenland ice sheet.  Qualitatively, the Stommel box model (<ref>) mirrors certain aspects of the non-autonomous Stommel box model (<ref>) close to the NSF (by setting $A=0$ in the non-autonomous case). The former presents a relatively generic and simple model with a region of bistability of two stable states that lose stability via a smooth saddle-node bifurcation while the latter undergoes an NSF bifurcation. The non-autonomous Stommel box model provides a framework in which we can develop critical (tipping) points in three cases: slowly varying bifurcation parameter only, high frequency oscillatory forcing only, and both types of variation combined.

Using the setup, we tested whether it would be possible to reconstruct the true value of the model parameters by assimilating ocean temperature and salinity observation and use those parameters to forecast the future circulation regime of the AMOC. Our twin experiments showed that the DA method is indeed capable of reconstructing temperature diffusivity and the advection parameters from observations. It also bring the salinity diffusivity closer to the truth, but the 18 year DA period used is too short to convergence to the minimal obtainable error. Using observed temperature and salinity from the EN4 dataset, the most-likely values of the model parameters have been determined with the DA system indicating that uncertainties in these estimates are smaller than 16%. In the absence of climate change future temperature and salinity is expected to lie in the same range of temperatures and salinities that have already been observed in the period 2004-2022. With the effects of climate, salinity differences between equatorial and pole box can move outside the current range because of the freshening caused by the melting Greenland ice sheet. Running the model with different climate change scenarios showed that this freshening is relatively more important to the direction of the AMOC circulation than changes in surface temperature. In particular, with the currently observed ice sheet melting rates, the model does not expect a reversal of the circulation direction to occur.

Though an interesting proof-of-concept that data assimilation can help with reconstructing Stommel model parameters, the forecasts produced in this study should be interpreted with some care. The simplicity of the Stommel box model means that it cannot represent the wide range of processes occurring in the North-Atlantic, nor can it represent the North-Atlantic geometry. In particular, we have seen that our model underestimates the transports. This is possibly because the southward boundary of the model was chosen to lie too far north of the equator. Next to this, the two box model considers the North-Atlantic to be a closed system and does not account for the southward transport of North-Atlantic Deep Water past the equator towards the Antarctic Circumpolar Current nor for the northward flow of Antarctic Bottom Water and waters from the Angulhas current <cit.>. Accommodating these concepts requires more complex conceptual models like the Rooth model <cit.> or full climate models.

As a last note, the simplicity of the chosen reduced order model makes our analysis possible as we are able to initialize large ensembles at low computational costs. This setup can serve as a baseline in complexity with more detailed versions of the model available for further studies on higher order dynamics.

§ DATA AVAILABILITY

The MetOffice EN4 objective analysis temperature and salinity data with XBT and MBT corrections used to create observations for assimilation as well as the surface forcing is available at <https://hadleyserver.metoffice.gov.uk/en4/download-en4-2-2.html> as the EN.4.2.2.analyses.g10.XXXX.zip files. The hadley_obs.py script used to process these observations can be found in the inversion branch of the <cit.> GIT repository. Also included in this fork of DAPPER <cit.> are the code for the Stommel model (dapper/mods/stommel/__init__.py) as well as the scripts used to carry out the data assimilation (dapper/mods/stommel/ref_da_forcing.py) using these observationsand the scripts for the twin experiments.

§ ACKNOWLEDGEMENT

We would like to thank John Gemmer and Christopher K.R.T. Jones for starting and funding this research project during the 2023 MCRN Summer School and Research Program. Ivo Pasmans would like to acknowledge the support of the Scale-Aware Sea Ice Project (SASIP) funded by Schmidt Sciences – a philanthropic initiative that seeks to improve societal outcomes through the development of emerging science and technologies. Emmanuel Fleurantin was supported by NSF grant DMS-2137947. A part of this project originated from the Mason Experimental Geometry Lab (MEGL) at George Mason University. A special thanks goes out to those MEGL students who started with the twin experiment setup in the Fall 2022 semester.