A functional regression model for heterogeneous BioGeoChemical Argo data in the Southern Ocean

The rapid growth of high-resolution data motivates continued advances in flexible and efficient statistical methodology. One prime example is the Argo data, which is collected by autonomous devices drifting through the world oceans called floats that measure temperature and salinity. Argo floats collect their data in so-called profiles, collections of measurements at various depths for fixed points in time and space (see Figure 1 and Argo, 2020). The Argo data has greatly enhanced our understanding of the physical nature of the upper 2000 meters of the oceans with near-uniform spatial sampling of these variables (Johnson et al., 2022). Argo floats determine their own depth through pressure measurements as an increase in 1 decibar in pressure occurs for each additional meter of depth (approximately). From here on, we will thus often use pressure rather than depth to refer to a float’s vertical position in the ocean.

An important expansion of the Argo Project’s goals is Biogeochemical (BGC) Argo, which has added a limited number of floats that also measure oxygen and nitrate, among other variables. We focus on Argo data in the Southern Ocean because the majority of the BGC Argo data has been collected under the Southern Ocean Carbon and Climate Observations and Modeling (SOCCOM) project (Riser et al., 2023, soccom.princeton.edu). The increased availability of biogeochemical data is important to improve our understanding of vital biogeochemical processes such as the biological carbon pump and air-sea CO₂ exchanges, monitor changes such as ocean deoxygenation and acidification, and improve estimates of the carbon budget (Bittig et al., 2019; Claustre et al., 2020; Gray et al., 2018). However, BGC floats use additional sensors, leading to higher fixed and maintenance costs than Core Argo floats (that is, those that measure only temperature and salinity). Thus, it is not practical to introduce the same coverage as achieved by Core Argo floats. Indeed, in the Southern Ocean during 2020, Core Argo floats collected approximately 12000 profiles while BGC Argo floats collected just 3600 over the same time period (see Figure 1).

When modelling biogeochemical data, it is thus desirable to utilize the more prevalently available Core Argo data in addition to the available BGC data to improve prediction and estimates and reduce uncertainty. This paper is aimed towards this goal, and our main contributions are as follows.

1. 1.

   We develop a descriptive statistical model on the joint distribution of oxygen, temperature, and salinity in the Southern Ocean. The model fully takes into account the functional nature and space-time information of the Argo data.

2. 2.

   Based on this model, and leveraging the more readily available temperature and salinity data, we address methodological and computational issues related to the functional prediction and interpolation of oxygen levels at arbitrary locations.

The basic building block for our model is a spatio-temporal, function-on-function regression model that encodes the functional relationship between the temperature, salinity and oxygen profiles across all pressure levels ranging from $0$ to $2000$ decibars. This enables us to predict oxygen profiles using temperature and salinity profiles. It also facilitates the estimation and prediction of other quantities of scientific interests such as derivatives and integrals of the profiles (Park et al., 2023).

One major challenge of modelling Argo data in the Southern Ocean is that the transition from warmer subtropical waters in the north to colder Antarctic waters does not occur smoothly. Instead, ocean properties in that region are characterized by a series of sudden changes, occurring at so-called fronts — sharp transition zones generally aligned east–west (Chapman et al., 2020; Orsi et al., 1995). In particular, the contrasting behavior of oxygen in various regimes of the Southern Ocean has been well established in Bushinsky et al. (2017). Consequently, any stationarity assumption on the data across different fronts would be violated. However, it is established that ocean properties in any region between fronts tend to be homogeneous (Chapman et al., 2020). To address this challenge, we augment our functional regression model with a latent mixture component aimed at capturing regional information. Our model thus describes aspects of the joint distribution of oxygen, temperature, and salinity modulated by different regimes in the Southern Ocean. In our analyses, we demonstrate that the dependence structure in pressure of these three variables varies critically across the regimes estimated through the mixture model. As a byproduct we obtain estimates for the locations of ocean fronts and how they vary with time. These are of considerable scientific interest in their own right and have been the subject of research in the oceanographic literature (Chapman et al., 2020). Consequently, they constitute an additional key contribution of this work.

Identifying ocean fronts through the use of clustering and mixture models in oceanography has recently gained significant traction. For example, Jones et al. (2023, 2019) focus on the Southern Ocean region; Rosso et al. (2020), Boehme and Rosso (2021), and Fonvieille et al. (2023) consider subregions of the Southern Ocean, while Maze et al. (2017), Poropat et al. (2023), and Houghton and Wilson (2020) analyze other oceanographic areas. Most of these works use standard Gaussian mixture models. The primary difference between this literature and our work here is that we account for the spatio-temporal dependence for the cluster memberships and the data within each cluster while these critical aspects are ignored in the above cited works. We accomplish this by encoding the spatio-temporal dependence of cluster memberships using a Markov random field model (Liang et al., 2020; Jiang and Serban, 2012; Besag, 1975). In particular, this enables predictions at any location whether temperature and salinity data are available or not. An important application of this is mapping, or the creation of gridded products, where data sampled at irregular locations in space and time are used to provide estimates at a standard grid (for example, Roemmich and Gilson, 2009; Gray and Riser, 2014).

When temperature and salinity measurements are available, predictions of oxygen may be straightforwardly made using off-the-shelf methodology. A well-known example in that regard is the random forest approach of Giglio et al. (2018), which produces highly accurate predictions of oxygen at fixed pressure levels using temperature and salinity data. For this scenario, we provide comparisons of our mode-based approach with the random forest approach of Giglio et al. (2018). However, when estimating the oxygen level at a location or a pressure level where temperature and/or salinity information is not available, the random forest approach requires an additional step of interpolating predicted oxygen data onto this location. Our model-based approach, on the other hand, readily provides direct predictions as well as estimated uncertainty of oxygen levels at arbitrary locations and pressure levels using nearby temperature and salinity data without interpolations.

Our approach also involves advancements in both statistical methodology and implementation. There is a rich literature on function-on-function regression, including the seminal paper Yao et al. (2005) and Zhou et al. (2008), while modelling spatio-temporal dependence in univariate functional data has been studied extensively by, e.g., Zhou et al. (2010); Jiang and Serban (2012); Liang et al. (2020). However, we are unaware of any development of function-on-function regression models that accounts for spatio-temporal correlation in the functional response and predictor variables. Furthermore, to the best of our knowledge, functional kriging and cokriging problems in the context of functional mixture models have not been explored at all. The SOCCOM Argo data considered here alone consists of more than fifteen thousand locations and millions of observations. Consequently, reducing the runtime and memory usage of the model estimation and prediction procedure is essential. This work is accompanied by an R package implementing our methodology where crucial parts of the estimation and prediction procedure are written in C++.

We finally briefly describe the organization of the rest of the paper. In Section 2, we propose a novel mixture model for multivariate, functional data with spatial dependence. In Section 3, we discuss parameter estimation, prediction, and model selection. In Section 4, we describe our analysis on the Argo data, including the description of the estimated model, evaluation of oceanographic front estimates, and cross-validation experiments. Finally, we discuss limitations and potential improvements to our work in Section 5.

Argo data is publicly available from the Argo Global Data Assembly (GDAC, Argo, 2020) after it undergoes various data control measures (Maurer et al., 2021). A precise description of the data under consideration is given in the supplement Section S1.

Notationally, our data consists of profiles at $n$ locations $\{x_{i}\}_{i=1}^{n}$ in the Southern Ocean and corresponding points in time $\{t_{i}\}_{i=1}^{n}$. We use $\mathsf{D}$ to denote the spatial domain and use the positive real axis $\mathbb{R}_{+}$ as time domain. To simplify notation, we will denote a space-time coordinate $(x_{i},t_{i})\in\mathsf{D}\times\mathbb{R}_{+}$ by $s_{i}$. At each point $s_{i}$, we observe profiles of a biogeochemical quantity $Y$ (oxygen) as well as temperature $T$ and salinity $S$ each at pressure levels $\{p_{ij}\}_{j=1}^{n}$. Consequently, oxygen data for one profile for $Y$ is then $\{Y_{s_{i}}(p_{ij})\}_{j=1}^{n_{i}}$ and for temperature and salinity $\{T_{s_{i}}(p_{ij})\}_{j=1}^{n_{i}}$ and $\{S_{s_{i}}(p_{ij})\}_{j=1}^{n_{i}}$. Even though the number of observations as well as the pressure indexes may vary between temperature, salinity and oxygen, we suppress this in our notation for clarity of presentation. Motivated by this structure in the data, we model the observations as point evaluations from unknown random functions of pressure, observed with additive measurement errors:

$\displaystyle Y_{s_{i}}(p_{ij})=\mathcal{Y}_{s_{i}}(p_{ij})+\epsilon_{ij}^{Y},\quad T_{s_{i}}(p_{ij})=\mathcal{T}_{s_{i}}(p_{ij})+\epsilon_{ij}^{T},\quad S_{s_{i}}(p_{ij})=\mathcal{S}_{s_{i}}(p_{ij})+\epsilon_{ij}^{S}.$

We use calligraphic letters to denote the underlying functions and regular letters for the observations. For each index $s$, we model the functions $\mathcal{Y}_{s},\mathcal{T}_{s}$ and $\mathcal{S}_{s}$ as elements of a Hilbert space $\mathbb{H}$ which we choose to be $\mathbb{H}=L_{2}([0,2000])$, the Hilbert space of square integrable functions on the closed interval $[0,2000]$ because most Argo profiles consist of measurements between 0 and 2000 decibars. In addition, we assume sufficient regularity of the functions to ensure pointwise evaluations, and we later impose a smoothness penalty on the 2nd derivative. The measurement errors $\{\epsilon_{ij}^{Y}\},\{\epsilon_{ij}^{T}\}$, and $\{\epsilon_{ij}^{S}\}$ are assumed to be i.i.d. Gaussian random variables with zero mean and variances $\sigma_{Y}^{2},\sigma_{T}^{2}$ and $\sigma_{S}^{2}$ to be estimated from the data. The main focus of our model is the biogeochemical quantity $Y$, so we collect temperature and salinity in a covariate vector of functions $\mathcal{X}_{s}(\cdot,\cdot)=(\mathcal{T}_{s}(\cdot),\mathcal{S}_{s}(\cdot))\in\mathbb{H}\times\mathbb{H}$. Combined, this yields two space-time, function-valued random fields $\{\mathcal{Y}_{s}:s\in\mathsf{D}\times\mathbb{R}_{+}\}$ and $\{\mathcal{X}_{s}:s\in\mathsf{D}\times\mathbb{R}_{+}\}$.

Modelling Ocean Fronts via a Latent Mixture Component As stated in the introduction, one main challenge with the Argo data is its nonstationary. To account for this, much of the previous work on the Argo data uses local approaches for example via sliding windows; see for example, Kuusela and Stein (2018); Yarger et al. (2022); Park et al. (2023). A locally-stationary model is estimated for each window, and since parameters are estimated anew for each window, the model adapts to varying covariance structures in different locations. However, in the Southern Ocean, this approach may not be suitable due to the sharp and sudden changes in water properties at the ocean fronts. Moreover, due to the circumpolar nature of the ocean fronts, water properties may stay homogeneous across long distances which any local approach will fail to fully utilize. In order to explicitly model the ocean fronts, we partition the data into $G$ clusters, with ocean fronts representing the boundaries between clusters. We use a latent field $\{Z(s)\}$ where each $Z(s)$ is a random variable taking value in $\{1,\dots,G\}$ for some $G$ to be specified. We then model

$\displaystyle\mathcal{Y}_{s}(p)=\sum_{g=1}^{G}\mathbb{I}(Z(s)=g)\mathcal{Y}_{s}^{g}(p),\quad\mathcal{X}_{s}(p_{1},p_{2})=\sum_{g=1}^{G}\mathbb{I}(Z(s)=g)\mathcal{X}_{s}^{g}(p_{1},p_{2})$

where the random fields $\{\mathcal{Y}_{s}^{g}\}$ and $\{\mathcal{X}_{s}^{g}\}$ govern the spatio-temporal covariance structure in cluster $g$. Here $\mathbb{I}$ is the usual indicator function.

For $\{Z(s)\}$, we employ a commonly used Markov random field model (cf. Jiang and Serban, 2012; Liang et al., 2020), where the marginal distribution $Z(s_{i})$ given appropriately defined neighbors $\partial s_{i}$ is given by

$\displaystyle\mathbb{P}\left(Z(s_{i})=g|\partial s_{i}\right)\propto\exp\bigg{(}\xi\sum_{s_{j}\in\partial s_{i}}\frac{\mathbb{I}(Z(s_{j})=g)}{\mathrm{dist}(s_{i},s_{j})}\bigg{)}.$

(1)

We use $k$-nearest neighbor graphs to define the neighborhood structure based on the distance $\mathrm{dist}(s_{i},s_{j})=\mathrm{dist}_{\mathrm{long}}+3\mathrm{dist}_{\mathrm{lat}}+12\mathrm{dist}_{\mathrm{dayofyear}}$ to account for the anisotropy in the Argo data. The additional weighting by the inverse distance follows the intuition that profiles that are close should provide better cluster information than profiles that are further away. This also reduces the importance of the choice of the number of nearest neighbors to construct the underlying graph which we chose to be 15. Consequently, we found that reasonably varying this number has negigible impact. The Markov random field model allows cluster memberships to be interpolated between the irregularly-sampled profile locations. Moreover the Markov structure yields computational advantages by using a pseudo-likelihood approximation as in Besag (1975).

We found that incorporating a temporal component in the distance is important. Indeed, we implemented a Markov random field that used only space (without a seasonal component) and found that in many locations this gave conflicting information on cluster membership from winter versus summer profiles.

Spatio-Temporal Function on Function Regression Within each cluster, we use the following spatial regression model to specify the joint distribution of $\mathcal{X}^{g}$ and $\mathcal{Y}^{g}$:

$\displaystyle\mathcal{Y}_{s}^{g}(\cdot)=\mu^{g}_{\mathcal{Y},t}(\cdot)+\mathscr{F}^{g}(\mathcal{X}_{s}^{g}-\mu^{g}_{\mathcal{X},t})(\cdot)+\mathcal{E}_{s}^{g}(\cdot),$

(2)

where the $\mathscr{F}^{g}:\mathbb{H}\times\mathbb{H}\to\mathbb{H}$ are linear operators and the $\mathcal{E}^{g}_{s}$ are functional, spatio-temporally correlated residuals that are independent of $\{\mathcal{X}^{g}_{s}\}$. Here, $\mu^{g}_{\mathcal{Y},t}(\cdot)$ and $\mu^{g}_{\mathcal{X},t}(\cdot)$ are mean functions of pressure for the response and predictors, respectively, in the $g$-th group and depend on time $t$. In the mean function $\mu^{g}_{\mathcal{Y},t}(\cdot)$, we include seasonality using a Fourier basis as in Roemmich and Gilson (2009):

$\displaystyle\mu^{g}_{\mathcal{Y},t}(\cdot)$

$\displaystyle=\gamma^{g}_{\mathcal{Y},0}(\cdot)+\sum_{l=1}^{R}\sin\left(2\pi l\frac{t}{365.25}\right)\gamma^{g}_{\mathcal{Y},l}(\cdot)+\sum_{l=1}^{R}\cos\left(2\pi l\frac{t}{365.25}\right)\gamma^{g}_{\mathcal{Y},R+l}(\cdot).$

The evaluation $\gamma^{g}_{\mathcal{Y},0}(p)$ then represents the average (across the year) mean of oxygen for group $g$ at pressure $p$, while the $\gamma^{g}_{\mathcal{Y},l}(p)$ for $l=1,\dots,2R$ describe seasonal variations in that mean based on sines and cosines. We take $R=3$ in this case to capture the main seasonality component, which is half of the number basis functions used in Roemmich and Gilson (2009) because in Antarctica there are only two seasons (NASA, 2024). The mean functions of the predictors have identical form. Next, we employ a Karhunen-Loève type expansion along the pressure dimension,

$\displaystyle\mathcal{X}_{s}^{g}(\cdot,\cdot)=\mu^{g}_{\mathcal{X},t}(\cdot,\cdot)+\sum_{l=1}^{\infty}\alpha^{g}_{l}(s)\psi^{g}_{l}(\cdot,\cdot),\quad\mathcal{E}^{g}_{s}(\cdot)=\sum_{l=1}^{\infty}\eta^{g}_{l}(s)\phi^{g}_{l}(\cdot)$

(3)

in order to model the spatio-temporal correlation structure through the principal component random fields $\{\alpha^{g}_{l}(s):s\in\mathsf{D}\times\mathbb{R}_{+}\}$ and $\{\eta^{g}_{l}(s):s\in\mathsf{D}\times\mathbb{R}_{+}\}$ (Zhou et al., 2010). The functions $\{\psi^{g}_{l}\}_{l=1}^{\infty}$ and $\{\phi^{g}_{l}\}_{l=1}^{\infty}$ form an orthonormal basis for $\mathbb{H}\times\mathbb{H}$ and $\mathbb{H}$ respectively for each $g$ and are called the principal component functions (PCFs). In practice, the infinite sums have to be truncated appropriately. Through the joint principal component decomposition of $\mathcal{X}_{s}^{g}(\cdot,\cdot)$, dependence between predictors (temperature and salinity) is modeled. Quantifying this dependence in the model is critical, since these variables are known to be dependent. For example, temperature, salinity, and pressure jointly determine potential density which is known to be generally monotone as a function of pressure (Section 3.5 of Talley, 2011). The coefficient random fields $\{\alpha^{g}_{l}(s):s\in\mathsf{D}\times\mathbb{R}_{+}\}$ and $\{\eta^{g}_{l}(s):s\in\mathsf{D}\times\mathbb{R}_{+}\}$ are assumed to be independent across cluster $g$ and component $l$. Using (3) in (2), and afterwards truncating the infinite sums after $Q_{1}$ and $Q_{2}$ terms, we obtain

$\displaystyle\mathcal{Y}^{g}_{s}(\cdot)=\mu^{g}_{\mathcal{Y},t}(\cdot)+\sum_{l=1}^{Q_{1}}\alpha^{g}_{l}(s)(\mathscr{F}^{g}\psi^{g}_{l})(\cdot)+\sum_{l=1}^{Q_{2}}\eta^{g}_{l}(s)\phi^{g}_{l}(\cdot).$

We model the random fields constituting the principal component scores as Gaussian random fields as is common in spatial statistics. Because they have zero mean, it remains to specify their covariance structure. A popular choice is the Matérn covariance function as advocated for by Stein (2013). However, for the Argo data an isotropic covariance function such as the standard Matérn covariance is likely not appropriate. For example, Kuusela and Stein (2018) found that, typically, the correlation range scales are different for longitude and latitude. We follow Guinness (2021) who, when modelling Argo data, used spatial deformation of the sphere via gradients of the spherical harmonics to model anisotropies. Concretely, we deform the sphere by mapping a point $x\mapsto x+\Phi_{g,\alpha,l}(x)$ where $\Phi_{g,\alpha,l}(x)$ is the weighted sum of gradients of the spherical harmonics of degree 2 at $x$. The weights will be learned from the data. We then use the following covariance

$\displaystyle\mathbb{C}\textrm{ov}\{\alpha_{g,l}(s),\alpha_{g,l}(s^{\prime})\}=\sigma^{2}_{g,\alpha,l}\left\{d_{g,\alpha,l}(s,s^{\prime})\right\}^{\nu^{g,\alpha,l}}K_{\nu^{g,\alpha,l}}\{d_{g,\alpha,l}(s,s^{\prime})\},$

(4)

where $\nu\geq 0$, $K_{\nu}$ is the modified Bessel function of the second type and for $s=(x,t)$,

$\displaystyle d_{g,\alpha,l}(s,s^{\prime})=\frac{\lVert x+\Phi_{g,\alpha,l}(x)-\{x^{\prime}+\Phi_{g,\alpha,l}(x^{\prime})\}\rVert}{\kappa_{g,\alpha,l,x}}+\frac{\lvert t-t^{\prime}\rvert}{\kappa_{g,\alpha,l,t}},$

where $\|\cdot\|$ is the standard Euclidean norm (chordal distance). An R implementation of this covariance function is readily available in the GpGp package (Guinness et al., 2021). Generally, we found that allowing more flexibility in the covariance function improved prediction performance. The smoothness parameter $\nu^{g,\alpha,l}$, variance parameter $\sigma_{g,\alpha,l}^{2}$, and range parameter $\kappa_{g,\alpha,l,x}$ and $\kappa_{g,\alpha,l,t}$ govern properties of the Matérn model for the $l$-th predictor principal component of the $g$-th group. We use the same covariance form for $\{\eta_{g,l}(s)\}$.

The mean functions, the PCFs, and the functions $\{\mathscr{F}^{g}\psi^{g}_{l}\}$ are unknown and have to be estimated from the data. To do so, we use a cubic B-spline basis, a common choice in functional data analysis with favourable computational and provable optimality properties (Wahba, 1990; Hsing and Eubank, 2015). Specifically, if $\{b_{i}\}_{i=1}^{P}$ are the elements of the cubic B-spline basis on [0, 2000] write $B(p)=(b_{1}(p),\dots,b_{P}(p))$. Then we use

	$\displaystyle\gamma_{\mathcal{Y},l}^{g}(p)=B(p)\Upsilon_{\mathcal{Y},l}^{g},\quad\phi^{g}_{l}(p)=B(p)\Theta^{g}_{\mathcal{E},l},\quad(\mathscr{F}^{g}\psi^{g}_{l})(p)=B(p)\Lambda^{g}_{l},$
	$\displaystyle\gamma_{\mathcal{X},l}^{g}(p_{1},p_{2})=(B(p_{1}),B(p_{2}))\Upsilon_{\mathcal{X},l}^{g},\quad\psi^{g}_{l}(p_{1},p_{2})=(B(p_{1}),B(p_{2}))\Theta^{g}_{\mathcal{X},l},$

where $\Upsilon_{\mathcal{Y},l}^{g}$ is the $l$-th column of the $P\times 7$ matrix of spline coefficients for the seasonal means of $\mathcal{Y}_{g}$ and $\Upsilon_{\mathcal{X},l}^{g}$ is the corresponding matrix for the covariates, $\Theta^{g}_{\mathcal{E},l}$ is the $l$-th column of the $P\times Q_{2}$ matrix containing the coefficients for the PCFs for $\mathcal{E}$, and the $P\times Q_{1}$ matrix $\Lambda^{g}$ contains the coefficients for the $\{\mathscr{F}^{g}\psi^{g}_{l}\}$. Similarly, $\Theta^{g}_{\mathcal{X}}$ is a $2P\times Q_{1}$ matrix. For the PCFs, we impose the following constraints to enforce orthonormality:

$\displaystyle\left(\Theta_{\mathcal{X}}^{g}\right)^{\top}\left(I_{2}\otimes\int B(t)^{\top}B(t)\,\mathrm{d}t\,\right)\Theta_{\mathcal{X}}^{g}=I_{Q_{1}},~{}~{}~{}~{}\left(\Theta_{\mathcal{E}}^{g}\right)^{\top}\int B(t)^{\top}B(t)\,\mathrm{d}t\,\Theta_{\mathcal{E}}^{g}=I_{Q_{2}}.$

(5)

Here $\otimes$ is the standard Kronecker product and $I_{q}$ is the $q\times q$ identity matrix. For identifiability reasons, we require the first entries in each column of $\Theta_{\mathcal{X}}^{g}$ and $\Theta^{g}_{\mathcal{E}}$ to be positive and order the PCFs by the marginal variances of the random fields constituting the principal component scores.

Here, we present results from our analysis of the Argo data, beginning with the estimated model in Section 4.1. We then provide a cross-validation study to assess prediction performance in Section 4.2 and present our data product of comprehensive functional predictions in Section 4.3.

We have developed and implemented a comprehensive space-time functional regression mixture model for temperature, salinity, and oxygen. The model yields a functional cokriging framework, where the oxygen concentration as a function of depth and the season can be estimated over the entire Southern Ocean along with (functional) uncertainties. The model has been validated against the most accurate existing non-functional approaches, and the mixture component of the model identifies contrasting distributions of the measured variables.

While our analysis introduces a comprehensive study of the dependence structure and prediction of oxygen in the Southern Ocean, there are a number of areas where the richness and complexity of the Argo data calls for improvements. As the principal components are entirely based on the pressure domain, predictions may be better if the components considered spatial information as well (cf., Kuenzer et al., 2020). A nonfunctional prediction approach, as done in Salvana and Jun (2022), might also fare well, but requires more computing than readily available at the scale of the entire Southern Ocean. Our model also does not take into account the monotonicity of density or the soft constraint on the saturation of oxygen near the ocean surface, which could improve predictions; however, applying standard functional data analysis approaches were not successful due to the limited domain in the density dimension for some profiles. Addressing this problem with the recent research of (for example) Lin and Wang (2022) may be possible. The current model is also developed with independence assumed between different clusters; incorporating dependence between clusters may better represent the available data. This, as well as evaluating whether independence should be assumed for different principal components the principal component scores may be assumed independent (as in Liang et al., 2023) would more comprehensively justify the data analysis. Including scalar covariates like satellite data in the model is straightforward and may be of practical interest. Finally, using more sophisticated model selection for the Markov random field might improve the cluster predictions.

Statistical approaches like the one developed here will play an important role for the analysis of data from the planned global deployment of BGC Argo floats (see NSF Award 1946578).

All data used in this work is publicly available. The Argo data is available at Argo (2020), and we use the SOCCOM-specific data files from Riser et al. (2023). For plotting of traditional estimates of fronts, we use the sea ice data Fetterer et al. (2017) and the Roemmich-Gilson mean Roemmich and Gilson (2009). Finally, we use the grid for prediction from the Southern Ocean State Estimate (Verdy and Mazloff, 2017). Code to produce results presented in the paper are available at https://github.com/dyarger/argo_fstmr_oxygen/.

Data were collected and made freely available by the Southern Ocean Carbon and Climate Observations and Modeling (SOCCOM) Project funded by the National Science Foundation, Division of Polar Programs (NSF PLR-1425989 and OPP-1936222), supplemented by NASA, and by the International Argo Program and the NOAA programs that contribute to it. (https://argo.ucsd.edu, https://www.ocean-ops.org). The Argo Program is part of the Global Ocean Observing System. Computational resources for the SOSE (Verdy and Mazloff, 2017) were provided by NSF XSEDE resource grant OCE130007 and SOCCOM NSF award PLR-1425989. This research was supported in part through computational resources and services provided by Advanced Research Computing at the University of Michigan, Ann Arbor. Authors Korte-Stapff, Stoev, and Hsing are supported for this work by NSF Grant DMS-1916226. Author Yarger is supported for this work by NSF Grant DGE-1841052.