{'fraction_non_alphanumeric': 0.04508661029472129, 'fraction_numerical': 0.015121911173138495, 'mean_word_length': 4.406657789613848, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 7, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 32, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A multivariate functional-data mixture model for spatio-temporal data: inference and cokriging

In this paper, we introduce a model for multivariate, spatio-temporal functional data. Specifically, this work proposes a mixture model that is used to perform spatio-temporal prediction (cokriging) when both the response and the additional covariates are functional data. The estimation of such models in the context of expansive data poses many methodological and computational challenges. We propose a new Monte Carlo Expectation Maximization algorithm based on importance sampling to estimate model parameters. We validate our methodology using simulation studies and provide a comparison to previously proposed functional clustering methodologies. To tackle computational challenges, we describe a variety of advances that enable application on large spatio-temporal datasets. The methodology is applied on Argo oceanographic data in the Southern Ocean to predict oxygen concentration, which is critical to ocean biodiversity and reflects fundamental aspects of the carbon cycle. Our model and implementation comprehensively provide oxygen predictions and their estimated uncertainty as well as recover established oceanographic fronts.


Introduction
The rapid growth of high-resolution and multivariate data motivates continued advances in flexible and efficient statistical methodology. One prime example is the Argo data, which is collected by floats that measure temperature and salinity at different depths in the ocean (see Figure 1 and Argo, 2020). 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 (Claustre et al., 2020). Biogeochemical floats require additional sensors leading to higher fixed and maintenance costs than Core Argo floats (Bittig et al., 2019). Thus, it is not practical to introduce the same coverage as achieved by Core Argo floats; in the Southern Ocean during 2020, Core Argo floats collected approximately 12,000 profiles (each being a collection of temperature and salinity measurements at varying depths for one location) while BGC Argo floats collected just 2,600 over the same time period (see Figure 1). However, biogeochemical information plays a vital role in the understanding of the world's carbon cycle as well as ocean ecosystems (Claustre et al., 2020). Consequently, the imbalance between Core Argo and BGC Argo floats calls for statistical methodology exploiting correlation of biogeochemical data with temperature and salinity to reconstruct biogeochemical data for Core Argo floats. One approach that has been employed by oceanographers is to fit a random forest, which represents well the complex relationships between the variables (Giglio et al., 2018). However, their analysis is restricted to one fixed depth at a time and consequently cannot model correlation between measurements stemming from the same profile. Instead, we view the Argo profiles as functional data indexed by space and time (Yarger et al., 2022). This approach leverages the correlation in depth and provides principled prediction of entire functions, dimension reduction through functional principal components analysis for efficient computation, and a concrete model for the covariance structure of the data.
One major challenge of modeling oceanographic 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 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). Consequently, any stationarity assumption on the data across different fronts will be violated. However, it is established that water properties between fronts tend to be homogeneous (Chapman et al., 2020). To handle the sudden changes in the distribution of the data caused by fronts while simultaneously exploiting the homogeneity of water properties between fronts, we propose a functional spatio-temporal regression mixture model.
Our work builds upon a strong foundation of research in functional data analysis and spatial statistics. To start, function-on-function regression has been comprehensively researched. We mention the seminal paper Yao et al. (2005) as well as Zhou et al. (2008) who formulate a likelihood-based estimation strategy for paired functional data. Modeling and estimation for spatially-correlated univariate functional data was studied in Zhou et al. (2010) and more recently in Zhang and Li (2021) who also considered piece-wise stationary functional data to accommodate spatial discontinuities in housing prices due to neighborhood characteristics. Spatial/spatio-temporal functional data that exhibit sudden changes are common, and substantial work has gone towards developing appropriate models. For example, Zheng and Serban (2018) investigate geographic disparities in pediatric health, Margaritella et al. (2021) model regions of brain activity, and Frévent et al. (2021) cluster abnormal unemployment rates. Recently, Liang et al. (2020) proposed a clustering methodology for univariate data allowing for cluster-dependent covariance structures. In a finite dimensional setting, Bolin et al. (2019) developed latent Gaussian random field mixture models for data observed on a regular grid, but the extension to a functional setting is not straightforward.
In the oceanographic literature, approaches for identifying fronts include statistical modeling, though most use Gaussian mixture models that ignore spatial information (Thomas et al., 2021;Chapman et al., 2020).
We describe next our contributions and how the paper is organized. In Section 2, we propose a novel mixture model for multivariate, functional data with spatial dependence.
We model the relationships between a univariate functional response and an arbitrary number of functional covariates. Furthermore, our approach can effectively handle irregularly-sampled data in both the spatial and functional dimensions of the data.
In Section 3, we discuss parameter estimation, prediction, and model selection. We introduce a Monte Carlo Expectation Maximization (MCEM) algorithm (Wei and Tanner, 1990) using importance sampling, which improves upon strategies employed in previous works (cf. Liang et al., 2020). In Section 3.3, we discuss how the estimated model can be used for predicting the response with available predictor data, while addressing the uncertainty of prediction based on that of the cluster memberships. Thus far, functional kriging and cokriging in the context of functional mixture models have not been studied. Furthermore, we provide a tractable AIC criterion in the presence of latent variables for model selection. In Section 4, we show in simulation studies that our EM algorithm produces viable parameter estimates and clustering results. In Section 5, we describe our analysis on the Argo data, which shows that the model provides quality predictions of biogeochemical variables with sensible uncertainty quantification. Furthermore, we recover established oceanographic fronts as a part of the clustering procedure. Finally, we discuss the outlook and potential improvements to our work in Section 6. The Argo data considered here consists of more than ten 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++, available at https://github.com/dyarger/argo_clustering_ cokriging. Our software seamlessly integrates a Vecchia approximation of the likelihood when modeling spatial dependence (Guinness, 2021). Furthermore, we exploit computational advantages of sparse matrices by circumventing the use of orthogonalized B-splines, an approach widely used in the functional data literature.
In this section we provide the framework for our modeling. Let S and T be two index sets.
The set S may correspond to a spatial index and the set T a temporal index, but this is not necessary. In our motivating application, S will correspond to a spatio-temporal index, and T will correspond to ocean depths. For a location s i ∈ S, where i = 1, 2, . . . , n, we observe measurements for a response Y and covariates X 1 , . . . , X K : The indices t ij and the number of observations n i may vary for different locations s i and across the covariates and the response. Our implementation reflects this, but for notational simplicity we omit their dependence on s i , X k , or Y . In our application, the response Y is a biogeochemical variable, while temperature and salinity are the covariates X 1 and X 2 .

Mixed Functional Spatio-temporal Linear Regression
We model sudden changes in the spatio-temporal distribution via an unobserved latent field where Z(s i ) takes values in {1, . . . , G} for some G to be identified. We will use the random field {Z(s)} to partition the data into G clusters and fit separate models for each cluster. To model spatio-temporal dependence, we assume that a Markov random field (MRF), where the distribution of Z i conditional on its neighbors ∂Z i is the Gibbs distribution with parameter ξ, that is This model was proposed in Jiang and Serban (2012) and used in Liang et al. (2020). We use graphs defined by k-nearest-neighbors ignoring the direction of the edges to obtain a symmetric neighborhood structure.
The observations are modeled as samples from unknown continuous random functions, observed with additive measurement errors. Notationally, we use calligraphic letters to denote the underlying functions and regular letters for the observations. Specifically, let where Y g = {Y g (s, t)} and X g,k = {X g,k (s, t)} are latent stochastic processes governing the spatio-temporal dependence within each cluster g. The cluster memberships Z i are unobserved. We assume that the underlying functional processes Y g and X g,k for different clusters are independent in g. The measurement errors for each variable are modeled as i.i.d. mean-zero Gaussian random variables.
For each cluster g, the {Y g (s i , ·)} and {X g,k (s i , ·)} are modeled as stationary functionvalued random fields taking values in the Hilbert spaces of functions H Y and H k respectively.
For simplicity, we assume that H Y = H 1 = · · · = H K =: H, and we choose H to be the space of square integrable functions on some compact interval [a, b]. We collect all of the covariates in a vector X g (s i , ·) := (X g,1 (s i , ·), . . . , X g,K (s i , ·)) , such that {X g (s, ·) : s ∈ S} is a random field taking values in H K , the direct sum of the Hilbert spaces. We endow H K with the weighted inner product f, g H K = K k=1 w k f k , g k H (Happ and Greven, 2018).
Here, µ g,Y and µ g,X are the mean functions of the random fields {Y g (s, ·)} and {X g (s, ·)}, respectively, F g is a bounded linear operator from H K to H, and the {E g (s, ·) : s ∈ S} form an H-valued, stationary, zero-mean random field representing a correlated spatiotemporal random effect in cluster g, which is independent from the covariates X g,k . The linear operator F g describes how the predictor processes relate to the response process, while E g represents variability unexplained by X g .

Dimension Reduction With Principal Components
Directly modeling functions in infinite-dimensional spaces is infeasible, so suitable representation and dimension reduction is necessary. To this end, we use a Karhunen-Loève expansion (cf. Hsing and Eubank, 2015) to write where µ X ,g and ψ g, take values in R K and φ g, takes values in R. denote the number of principal components for X g and E g respectively, our model becomes η g, (s)φ g, (t).

Modeling Spatio-temporal Correlation
Following Yarger et al. (2022), we model the spatio-temporal correlation of the random functions via the spatio-temporal correlations of their principal component scores. More precisely, we assume that the {α g, (s)} are stationary, zero-mean, Gaussian random fields that are uncorrelated across g and . Similarly, we assume that {η g, (s)} have the same structure and are independent of all {α g, (s)}. While the stationarity assumption seems restrictive, it is only imposed within each cluster g, and the combined model accommodates heterogeneity. For a specific data application, it remains to choose a suitable (valid) family of covariance function for the Gaussian processes. A commonly-used family is the Matérn covariance (Stein, 2013), which provides flexibility for the smoothness and range of spatial correlation in the process, written in closed form for parameters σ 2 , κ, ν > 0 as where K ν is the modified Bessel function of the second kind with order ν, and Γ is the gamma function. Here and throughout, · denotes the Euclidean norm on the respective space. The parameters are allowed to vary for different principal components and clusters.
We will describe the specific covariance functions we use in Sections 4 and 5.
3 Model Estimation

Function Approximation Using Polynomial Splines
To approximate the mean functions and the principal component functions, we use spline bases. For notational ease, we use the same spline basis for all functions, though this is not necessary. Let {b 1 , . . . , b P } denote the functions forming a P -dimensional cubic spline basis on the interval [a, b]. Let B(t) ∈ R 1×P denote the row vector of spline evaluations (b 1 (t), . . . , b P (t)). Then we approximate µ g,Y (t) ≈ B(t)θ g,Y where θ g,Y is a P -dimensional vector of spline coefficients. For the multivariate functional covariates, we write µ g,X (t) ≈ (I K ⊗ B(t)) θ g,X ,1 · · · θ g,X ,K =: where ⊗ is the standard Kronecker product and I K denotes the K × K identity matrix. Using similar notation, we have (ψ g,1 (t), . . . , ψ g,Q 1 (t)) ≈ (I K ⊗ B(t)) Θ g,X for the predictor principal component functions, ((F g ψ g,1 )(t), . . . , (F g ψ g,Q 1 )(t)) ≈ B(t)Λ g , and (φ g,1 (t), . . . , φ g,Q 2 (t)) ≈ B(t)Θ g,E , where Θ g,X , Θ g,E , and Λ g are matrices of spline coefficients of dimensions KP × Q 1 , P × Q 2 , and P × Q 1 . If α g i = (α g,1 (s i ), . . . , α g,Q 1 (s i )) and η g i = (η g,1 (s i ), . . . , η g,Q 2 (s i )) , the model for the data can be succinctly written as follows: Since the principal component functions are orthonormal, we impose the following orthonormality constraints for each g: Remark 1. A common approach to enforce the orthogonality constraints in (4) is to use an orthogonalized spline basisB(t) ∈ R 1×P and then require Θ g,X Θ g,X = I Q 1 and Θ g,E Θ g,E = I Q 2 ; see Zhou et al. (2010) and Liang et al. (2020), among others. However, for datasets as large as the Argo data, this becomes impractical because the evaluation matricesB(t) (at observed data points) for orthogonalized splines are less sparse compared to regular Bsplines. We outline in Section S5 how regular B-splines can be used to obtain algebraically identical estimates with great computational gains.
Remark 2. The model as defined in (2) and (3) is not identifiable as the eigenfunctions are unique only up to their signs. However, identifiability can be effectively enforced in estimation by requiring the first entries in each column of Θ g,X and Θ g,E to be positive, and ordering the eigenfunctions by the marginal variances of the principal component scores.

Monte Carlo EM Algorithm
We aim to fit the model parameters denoted Ω by penalized likelihood maximization.
First, we introduce some notation. We let Z := (Z 1 , . . . , Z n ) be the vector of cluster memberships, Y i ∈ R n i the response observations at s i , and X i ∈ R Kn i the vector of covariate observations arranged by covariate. Then we let Y := (Y 1 , . . . , Y n ) and X := (X 1 , . . . , X n ) be vectors of length n i=1 n i and K n i=1 n i , respectively. Define the random vector α i := G g=1 I(Z i = g)α g i . Finally, we use α g := (α g 1 , . . . , α g n ) ∈ R nQ 1 and α := (α 1 , . . . , α n ) ∈ R nQ 1 and define η, η g , and η i similarly.
To estimate Ω, we maximize the penalized log-likelihood of the observed data where P(Ω) denotes a roughness penalty and L denotes the likelihood function. We use the integrated squared second derivatives of the respective functions, which motivates our use of cubic B-splines (see Theorem 6.6.9 of Hsing and Eubank, 2015). For example, for the mean function for Y g and a smoothing parameter λ θ Y to be chosen, we use the penalty We use similar penalties for the principal component functions as well as for the functions {F g ψ g, }, and let P(Ω) then denote the sum of the individual penalties.
To compute Q Ω|Ω (m) for the E-step, using the law of total expectation, we have where Z denotes the set of all possible configurations of Z values. Consequently, this is a sum over G n terms which, even for small n, is computationally intractable since the likelihood of the complete data does not factor due to spatial correlation. Thus, we use a Monte Carlo EM algorithm (Wei and Tanner, 1990) as outlined next.
Monte Carlo E-Step. We propose to use importance sampling to estimate the expectation in Q Ω|Ω (m) . That is, we fix T ∈ N and create samples α (τ ) , η (τ ) , Z (τ ) for τ = 1, . . . , T from a proposal distribution h and approximate Q Ω|Ω (m) viâ where the quantities w the importance sampling weights. To simplify the notation, we will drop the dependence on Ω (m) below. For detailed information on importance sampling, see Robert and Casella (2004). Our goal is to construct an h that is easy to sample from and that is . First, note that for a given sample Z (τ ) , sampling from f (α, η|Z (τ ) , X, Y ) is tractable; see Proposition 1 in the supplement S1. It remains to find a proposal denoted P I for sampling Z (τ ) |X, Y and employ the following procedure: 1. Sample Z (τ ) from P I and compute w (τ ) .

Sample
where z = {z i } n i=1 and i ∼ j indicates that i and j are neighbors. Sampling from (8), known as the Potts model, has been studied in Feng et al. (2012) and references therein. When the vectors (X i , Y i ) are high-dimensional, the so-called external field term n i=1 log(f (X i , Y i |Z i = z i )) will dominate the other term of (8) in size. Due to this, we use a Gibbs sampler based on a single site updating scheme similar to Feng et al. (2012) to sample from P I (Z|X, Y ). With this, the final proposal distribution is h(α, η, Z) = P I (Z|X, Y )f (α, η|Z, X, Y ). The importance sampling weights are then using Bayes' Theorem, and the unknown constant C cancels in (7).
An alternative to importance sampling that has been considered is Gibbs sampling, see, for example, Liang et al. (2020). To describe this procedure, start with an initial sample of Z (0) , and then at the τ -th iteration generate samples via where sampling from P Z α (τ ) , η (τ ) , X, Y remains non-trivial. In their work, Liang et al. (2020) show that this can work well. However, in our application, we found that the resulting Markov chain tends to mix slowly, which we speculate is due to having a large number of observations per function. We illustrate this point in simulation studies, and provide a more thorough discussion in Section S3 in the supplement.

M-step.
We update our parameter estimates as Ω (m+1) = arg max ΩQ Ω Ω (m) . Since the log-likelihood of the complete data factors into separate terms in (6), updates for the parameters are done sequentially. A closed form solution exists for most of the parameters, though for the spatial parameters and ξ (as defined in Equation 1), an iterative approach is necessary. The updating formulas are given in Section S4, and the M-step also includes an orthogonalization step described in Section S5. For EM algorithms, initialization can play an important role. We provide a detailed initialization procedure in Section S8.

Prediction
A primary goal of the methodology is to predict at new locations, denoted byš = {š i }ň i=1 , with potentially available new predictor dataX. To accommodate new locations, we appropriately reconstruct the graph underlying the Markov Random Field using andš. In our application, we use the same number of nearest neighbors to do so. However, one may use more neighbors if there is a large number of new locations. Commonly, one is interested in predicting functions of Y(š i , ·). As an example, consider the evaluation so that prediction reduces to predicting functions A g of the latent principal component scores. Similarly, other functions of interest such as integrals and derivatives can be written in the same manner. Thus, in the following, let A be a function A : To quantify uncertainty, we use the conditional variance Var(A A A Y(š, ·)|W ).
Since the presence of unobserved variables Z andŽ = {Z(š i )}ň i=1 complicates the expressions of these moments, we again employ importance sampling for prediction. Similarly to the E-step, we generate samples Z (τ ) ,Ž (τ ) ∼ P I Z,Ž W for τ = 1, . . . , T with corresponding normalized weights w (1) , . . . , w (T ) , and use Monte Carlo integration via The right hand side of (10) averages the kriging error under the distribution of clusters, while (11) describes the predictive variance induced by uncertainty in cluster memberships.
If A is linear, expressions for A A A Y(š, ·) and the approximation of the conditional variance Var (A A A Y(š, ·)|W ) are obtained straightforwardly from multivariate Gaussian moments ofα g ,η g |W , Z,Ž via representations like (9). In the non-linear case, one may employ a parametric bootstrap approach. We provide explicit formulas in Section S9. Finally, inference for Z,Ž can be carried out based on the samples (Z,Ž) (1) , . . . , (Z,Ž) (T ) .

Model Selection
We propose to jointly choose the number of principal components and the smoothing parameters via an AIC criterion. The number of clusters G is assumed fixed and chosen using domain knowledge. Defining a model selection criterion for G and the Markov random field structure is highly nontrivial; see Stoehr et al. (2016) and references therein. Since we impose a smoothing penalty, the number of knots plays a minimal role in model fit, as long as there are sufficiently many to approximate relevant functions (Zhou et al., 2008).
Because the observed data likelihood is intractable, previous works (Liang et al., 2020;Huang et al., 2014) have used an approximation of likelihood-based criteria introduced as a last resort by Ibrahim et al. (2008), who suggested it may perform poorly when there are a large number of latent variables. Instead, we employ an importance sampling framework to approximate the marginal density f (X, Y ) directly at the cost of one E-step. Define corresponding to mixture model assuming i.i.d. observations. If we sample Z (τ ) from P I 2 , we may approximate the marginal density of X, Y via importance sampling as follows: where the last equality follows from Bayes' Theorem. Each of the terms on the right of (12) can be computed directly, with the exception of the normalizing constant of P Z (τ ) . This term, denoted by C 2 , is constant for all samples, smoothing parameters, and number of principal components. Thus, (12) allows us to estimate the AIC of a model M as where D(M) denotes the degrees of freedom of M. Then, for two models M 1 and M 2 , the term log(C 2 ) cancels in the difference AIC(M 1 ) − AIC(M 2 ). Here, we use P I 2 because we can sample from it and evaluate it, while it is not simple to evaluate P I . To determine the degrees of freedom of the model, we take into account the smoothing penalty imposed on functions as in Wei and Zhou (2010). This approach should not be used to select the number of clusters G, as C 2 would no longer be constant across models.
We choose different smoothing parameters for {µ g,Y }, {µ g,X }, {φ g }, {ψ g }, and {F ψ g }; for simplicity, we assume that these parameters do not vary over the cluster or principal component index. To avoid training many different models, we choose the penalty parameters anew in iterations of the EM algorithm based on our AIC criterion. While this incurs additional computation per EM iteration, it is significantly faster than training many separate models and performs well in our application and simulation studies.

Simulation
We test the performance of the estimation strategy in two simulation studies. The first study assesses clustering performance compared to other approaches. In the second study, we illustrate that our Monte Carlo EM algorithm can correctly estimate the model param- Cov (α g, (s), α g, (s )) = σ 2 g,α, exp (− s − s /κ g,α, ) with σ 2 g,α, and κ g,α, > 0. Both simulation studies are based on 100 independently-simulated datasets.

Clustering Accuracy
To assess our clustering strategy, we consider a setting comparable to Liang et al. (2020), who provide an extensive comparison of clustering methodologies for univariate functional data. We choose G = 2 groups, 3 principal components per group, and n = 200 locations.
The mean was taken as the same for both groups as µ(t) = e 2−2t cos(5(t − 0.1)), while the -th PC of the g-th group is taken as the (g + 2 )-th function of an orthogonal spline basis of dimension 11 with uniformly-spaced knots. We set the measurement error variance to be 1.0. We use a five-nearest-neighbors graph to describe the dependence structure of the Markov random field and set ξ = 0.5 in (1). The score variances were taken to be the same for both groups, with values of 1/3, 1/6, and 1/12, respectively. The spatial range parameters κ g, are defined by the matrix: κ = In Figure 2, we plot the percent of locations whose cluster memberships are correctly identified averaged over the simulations. In the dense setting (n i = 100), the Gibbs sampler from Liang et al. (2020) only slightly improves upon its initialization. This is because, across Gibbs samples and EM iterations, the sampled clusters do not change enough, indicating inadequate exploration of the set of possible cluster assignments. See Section S3 for more details. In contrast, importance sampling and the misspecified independent model appear to perform comparatively well in both scenarios with significant improvements over the initialization. This suggests that conditional on observations at a location s i , other (nearby) observations are only minimally informative when determining the cluster assignment Z i .
Thus, ignoring spatial correlation of the functional data in our proposal distribution appears

Parameter Estimation
We For the means, we use µ Y (t) = µ X 1 (t) = e 2−2t cos(5(t−0.1)) and µ X 2 (t) = 5 cos(5t), and we use an orthogonal cubic spline basis {φ 1 , . . . , φ 4 } based on two equidistant knots on [0, 1] as principal component functions for the spatio-temporal random effect. Finally, we and F ψ 4 (t) = 2e 1−t log(1−t+0.1)+1. The variances and spatial correlation parameters for the spatio-temporal random effect E are given in Table 1, along with their estimates. We select the penalty parameters and number of PCs based on (13), which chooses the correct number of response PCs 100% and predictor PCs 65% of the time. Figure 3 illustrates the resulting function estimates for µ Y , F ψ 1 , F ψ 4 , φ 1 , φ 2 , φ 4 . Overall, our MCEM algorithm appears to provide accurate estimates of the model parameters. Table 1 shows results for the variances and the spatial correlation scale parameters for the PCs of E. Zhang (2004) demonstrated that these parameters are not consistently estimable, however, consistent estimation of their ratio is possible and asymptotically, only the ratio matters for prediction.
Our estimates of the ratio are reasonably unbiased. Plots of the remaining functions and the corresponding table for σ α, and κ α, are in the supplement Section S6.

Data Analysis
In this section, we present the implementation and results on the Argo data in the Southern Ocean. We focus on this area since a large proportion of the biogeochemical ( (2020). We focus on oxygen data, though we also provide results for nitrate in Section S11.
Argo data quality control for oxygen is described in Maurer et al. (2021), and comprehensive analyses of oceanic oxygen include Keeling et al. (2010) and Stramma and Schmidtko (2021). For model training, we use about 13,000 BGC profiles which yields 1.6 × 10 6 oxygen measurements and 1.5 × 10 7 temperature and salinity measurements, while also building an auxiliary dataset of more than 150,000 Core profiles to be used for prediction in Section 5.3. In Section S7, we describe in detail how we created these datasets.
As discussed in the introduction, each Argo profile is a collection of measurements taken at different depths in the ocean. Thus, here the index t ij corresponds to pressure (depth) in the ocean, while s i represents the spatio-temporal index of a profile. The functional data perspective has some advantages in this context: integrals and derivatives of depth are immediately available, irregular sampling is handled by flexible B-spline bases, and dimension reduction in depth makes computation manageable. Viewing the Argo data as functional data has been considered in more detail in Yarger et al. (2022).
To provide background, we give a brief review of related oceanographic literature. First, interest in clustering applications in oceanography has risen quickly in recent years. Clustering has been used to characterize variability (Rosso et al., 2020), evaluate ocean heat content (Maze et al., 2017), examine El Niño events (Houghton and Wilson, 2020), and identify oceanographic fronts (Jones et al., 2019;Thomas et al., 2021;Chapman et al., 2020). Thus far, clustering methodology in oceanography has primarily consisted of Gaussian Mixture Models (for example, Rosso et al., 2020;Thomas et al., 2021) and k-means (Houghton and Wilson, 2020). These models have not taken advantage of spatial dependence of cluster memberships or observed data, which precludes spatial prediction.
As the BGC Argo program has developed, approaches for estimating and predicting its BGC variables has been an important area of research. Our aims closely match the spirit of Giglio et al. (2018), who use a random forest to predict oxygen from temperature, salinity, and spatio-temporal information. Expanding upon their work, we provide predictions at all pressure levels along with predicted uncertainties. In different settings, oceanographic approaches to predict BGC variables includes multivariate analysis (Liang et al., 2018), linear regression (Williams et al., 2016), and machine learning (Bittig et al., 2018).

Modeling and Implementation
Here, we specify our choice of spatial correlation function and other aspects of the application. For each principal component score, we use the anisotropic covariance function Cov(α g, (s), α g, (s )) = σ 2 g,α, exp − sin(γ(s, s )/2) 2 κ −2 g,α, ,s + |u − u | 2 κ −2 g,α, ,u , where γ(s, s ) is the angle in radians between warped locationss ands , and u and u are the respective times of the two samples. The warping is defined by gradients of spherical harmonic functions of degree 2 following Guinness (2021) and implemented in the R package GpGp . The parameters are the variance σ 2 g,α, > 0, the scales κ g,α, ,s and κ g,α, ,u > 0, and the parameters for the warping transformation. This is similar in spirit to the choice of correlation function in Kuusela and Stein (2018) and later in Yarger et al. (2022) because it allows the range of dependence to vary flexibly in space and time.
We chose 13 predictor and 11 response principal components using AIC. See results in Figure S6. We choose 5 clusters based on oceanographic evidence (Orsi et al., 1995;Chapman et al., 2020). For the Markov random field, we consider the 15 nearest neighbors based on the distance dist long + 3dist lat + 12dist day of year in kilometers and days, which takes into account seasonality and reflects scientific knowledge that scales in the longitude direction are often longer than latitude (for example, Roemmich and Gilson, 2009).

Hold-out Experiments
We assess the prediction performance of our model in two separate scenarios: leave-profileout and leave-float-out. For leave-profile-out, we hold out oxygen data from 1/5 of BGC profiles using a random partition, train the model on the remaining data, predict the left-out oxygen data, and repeat this process for of the five sets in the partition. For leave-float-out, we follow a similar process, but we instead partition floats into 10 sets, so that we hold out all oxygen data from selected floats at a time. In the leave-profile-out setting, nearby oxygen data from the same float is available; the leave-float-out setting is more challenging since one can only use oxygen data from other floats. Since Argo floats either have or do not have BGC sensors, leave-float-out is closer to the practical cokriging setting.
We found that for held-out locations with temperature and salinity data available, there is virtually no variation in cluster assignment across Monte Carlo samples used for prediction. Consequently, there are no mixture components, and it is straightforward leave-profiles-out (Right), summarized in 100 decibar bins. We compare predictions based on the cluster means only ("Mean"), the cluster means combined with temperature and salinity data ("Temperature and salinity"), and the full prediction ("Full prediction").
to construct confidence intervals using the resulting Gaussian distribution. We take the measurement error into account when constructing these confidence intervals.
In Figure 4, we plot the prediction error when using different components of the model.
Temperature and salinity are present at each left-out profile, so the predictions based solely on covariates are comparable between cross validation settings. The abundance of nearby oxygen data can only be used by the full model in the leave-profile-out setting. The resulting prediction error is close to the instrument measurement error around 3 µmol/kg (Maurer et al., 2021;Mignot et al., 2019). Figure 5 describes the performance of our uncertainty estimates. In the leave-profile-out setting, the intervals are a bit conservative; in the leave-float-out setting, they are slightly too narrow. Generally, however, the uncertainty estimates appear to work well. Thus, our model can adapt its uncertainty estimates to varying sampling settings, which is especially promising as BGC float coverage increases.
We compare our approach with the random forest approach of Giglio et al. (2018).
To assess prediction performance, we consider Argo measurements between 145 and 155 decibars and compare the point predictions for each held-out data experiment in Table 2.
For the random forest, we use the parameters described in Giglio et al. (2018) on these measurements. In the leave-profile-out setting, the random forest easily modulates nearby profiles into the predictions well, and both approaches are comparable to the instrument measurement error. In the leave-float-out setting, the functional approach provides comparable predictions. In addition, we note that out-of-the-box implementation of the random forest as used in Giglio et al. (2018) does not supply uncertainty estimates, and computational issues arose when we applied it to data across all pressure levels. Finally, our model allows inference on its components including the clustering structure and the functional linear model as discussed next in Section 5.3.  Table 2: Comparison with random forest (RF) and functional cokriging (FC, our approach here) at 145-155 decibars. The median absolute deviation (MAD) and root-mean-squarederror (RMSE) are the metrics used for each of the methods and leave-out settings.

Recovery of Oceanographic Fronts and Prediction Product
We now present results from our extensive analysis of Biogeochemical and Core Argo data in the Southern Ocean, which is based on interpolation of cluster memberships and oxygen on a grid using the methodology developed in Section 3.3. This allows us to provide interpretable maps of oceanographic fronts and oxygen concentration. First, we provide a thorough analysis of heterogeneity in ocean properties based on oxygen, temperature, and salinity data. In Figure 6, we plot predictions of cluster memberships and their uncertainty.
In addition, traditional estimates of fronts are plotted using the same criteria for fronts as Finally, we include the oxygen mean functions for each cluster.
We focus on cluster predictions for the austral summer (e.g. January) because historically, more measurements are available during this season, making the comparison between our results and existing literature more relevant. The boundary of our first cluster defines the Southern Antarctic Circumpolar front, which closely matches that found in Orsi et al.  Figure 6, while also characterizing using Equation (11) how this uncertainty contributes to that of the oxygen predictions, plotted in Figure S16.
To provide a set of oxygen predictions that can be used by scientists, we predict for each month of the year from July 2014 to June 2020. As an example in Figure 7, we plot oxygen predictions and their standard deviations at 150 decibars for January 15th, 2020, taking advantage of our functional approach that provides predictions at arbitrary pressure levels.
We then provide the predictions, their variances, and the averaged fields in a product in netcdf form, given in Section S10. We provide more prediction plots in Section S11.

Discussion
In this work, we introduce a model for multivariate and spatio-temporal functional data that enables large-scale cokriging (spatio-temporal prediction). The model simultaneously han-dles multivariate functional data, spatio-temporal dependence, and heterogeneous spatial features through a latent clustering structure. We introduce an importance-sampling-based Monte Carlo EM algorithm for model estimation that overcomes computational and inferential challenges faced by existing methodologies, as demonstrated in a simulation study. We contribute methodology for spatio-temporal prediction that takes into account uncertainty in cluster memberships. Among other methodological advances that make adaptation to large amounts of data possible, we avoid use of orthogonal B-splines, providing a framework that is more efficient in computation and memory. Our code and approach for estimation can be adapted for data in neuroscience or elsewhere in environmental sciences.
A number of extensions and applications of our methodology remain. For example, functional canonical correlation analysis may better represent the relationship between the covariates and the response. The current model is also developed with independence assumed between different clusters; incorporating dependence between clusters may better represent available data. 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.
We use our framework for a large-scale data analysis on the biogeochemical Argo data in the Southern Ocean. Via hold-out experiments and a comparison with existing approaches, we demonstrate that our model provides accurate predictions of oxygen data as well as reliable estimates of variability. Our clustering results independently recover established oceanographic fronts and provide insights into the spatio-temporal dynamics of these fronts.
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).

S1 Conditional Distributions
Here, we provide the conditional distributions required for our MCEM algorithm. The first proposition gives the conditional distribution of the observed data given the latent field Proposition 1. Let α be a normally distributed random vector of size m 1 with mean 0 and covariance Σ α . Let η be another normally distributed random vector of size m 2 with mean 0 and covariance Σ η . Furthermore, let X and Y by two random vectors with Y ∼ N (0, σ 2 Y I N 2 ) and X ∼ N 0, bdiag(σ 2 X 1 I n 1 , . . . , σ 2 X K I n K ) =: where the parameters σ 2 Y > 0 and σ 2 X k > 0 for k = 1, . . . , K are measurement error variances. Assume that α, η, X and Y are mutually independent of each other. Let N 1 = K k=1 n k and assume min(N 1 , N 2 ) > max(m 1 , m 2 ). Furthermore, let A 1 , A 2 , A 3 be three matrices of dimensions N 1 × m 1 , N 2 × m 1 and N 2 × m 2 respectively. Define Proof. The proof follows directly from standard properties of the normal distribution.
Next, we derive the conditional density f (α, η|X, Y , Z). As observations from different clusters are assumed to be independent, it suffices to derive the distribution for specific cluster. Thus to ease notation, we drop cluster dependence. For the response, let σ 2 Y = Var( Y 11 ). In addition, let Σ α = Var(α) and Σ η = Var(η) be defined by the spatial covariance.
Within the same cluster, we have Here, Θ E , Λ, and Θ X are the respective matrices for the principal component/random effect function coefficients and θ X and θ Y are the respective vectors for the mean coefficients, ignoring notation on the cluster.
Proof of Proposition 2. We have from which we obtain the inverse covariance matrix immediately. Let Σ inv α,α , Σ inv α,α and Σ inv η,η denote the blocks of Σ −1 . Then comparing the expression involving the means, we obtain and solving for µ α , µ η yields

S2 Computation of Conditional Distributions
Consider the notation in Section S1, and note that our goal is to sample from f (α, η|X, Y ) efficiently. All distributions should be considered conditional on Z, but this notation is omitted as in Section S1. Defining we may write In practice, Σ −1 will be considerably sparse due to the independence assumption on the scores; furthermore, one can apply popular sparse matrix spatial statistics approximations to A −1 , like Vecchia's approximation (Guinness, 2021). We next form the Cholesky matrix of Σ −1 = LL , which is the largest computation in the E-step. Then: 1. The conditional mean of α, η|X, Y can be found by backsolving and forwardsolving using L.

A mean-zero random normal vector with covariance Σ can be generated by backsolv-
ing against a vector with independent standard normal entries using L.
3. The marginal likelihood of X and Y (once again conditional on cluster memberships) can be computed, which factor into the importance sample weights. That is, write Finally, write n data = (K + 1) n i=1 n i as the number of total observed measurements, so that where |Σ data | denotes the determinant of Σ data . Using the matrix determinant lemma, The first determinant can be computed quickly from L, the second is quickly available when using a Vecchia approximation, and the final determinant (of a diagonal matrix) can also be computed quickly. Using the Woodbury matrix identity, The first term can be quickly computed, and the second term's only computational challenge is back/forward solving against L.
After tailoring this likelihood computation to a particular clustering Z, this solves the evaluation of f (X, Y |Z) in the importance sampling weights. Also, evaluation of f (X i , Y i |Z i ) can be done similarly.
These three computational derivations allow working with matrices the size of the total number of principal component scores, rather than the total number of data points. In practice, creating a permuted Cholesky matrix also improves sparsity and computation speed.

S3 Illustration with Gibbs Sampling
Here, we discuss the Gibbs sampling approach proposed by Liang et al. (2020) is poor if θ g and θ g , Λ g and Λ g , or Θ g and Θ g are substantially different. This is the case even if the covariance structures are the same in every cluster but the means are different.
and the respective quantities for the covariates, will be vanishingly small. This in turn will lead to high correlation of the samples of Z (τ ) ; that is, Z (τ ) may be equal to the initial value for all τ . This problem is particularly pronounced if the observation vectors X i and Y i are high-dimensional, as the conditional distribution of α will be increasingly concentrated around values representing the scores for cluster Z An alternative that avoids these issues is to use a Gibbs sampler by starting with an initial sample Z (0) and at the τ -th iteration proceed as follows: where now P Z α 1 (τ ) , η 1 (τ ) , . . . , α G (τ ) , η G (τ ) , X, Y is also a Potts model with external field. However, sampling f α 1 , η 1 , . . . , α G , η G Z (τ −1) , X, Y is computationally extremely expensive for large data sets, and in most cases P I is an appropriate proposal distribution, so that the importance sampling procedure provides a more efficient and comparably accurate approach to sample the clusters.

S4 Maximization Step Updates
In this section, we outline the parameter updates during the M-step. Throughout, we let = g , and P denote the relevant penalty matrices. Since we use selfnormalizing weights, we have T τ =1 w (τ ) = T , which are occasionally used interchangeably in the ensuing formulas. We denote the sampled predictor scores for the i-th profile from the τ -th Monte Carlo iteration as α , which correspond to principal components from the cluster Z i .
• The means for the response are updated via The updates for the means of the covariates are • We update the parameters in the matrices defining the principal components sequen-tially by column. The update for the -th column of Λ g is The update for the -th column of Θ X ,g is The update for the j-th column of Θ E,g is • The update for the measurement error of the response is Furthermore, for Θ X ,g and θ X ,g , let Θ k X ,g and θ k X ,g denote the rows from (k − 1)Q 1 + 1 through kQ 1 . Then, defining σ 2 X k = Var( k 11 ), we obtain an update via where X k,i are the measurements for profile i from the k-th predictor.
• Let α (τ ) g, be the random vector of observations of the random field of the -th principal component score for the covariates in cluster g from the τ -th sample from the importance sampling procedure. Let R (τ ) g, ,α denote the corresponding correlation matrix.
Furthermore, let n = g denote the number of sites belonging to cluster g for sample Z (τ ) . Correspondingly, let σ 2 g, ,α denote the marginal variance of the -th predictor score for the g-th group. Then the update for σ 2 g, ,α is The marginal variances of η for each group and principal component are updated in the same way.
• The number and kind of spatial (spatio-temporal) correlation parameters depends on the choice of correlation function. Regardless of that choice, we use numerical optimization of the relevant likelihood via the "L-BFGS-B" algorithm (Byrd et al., 1995) as implemented in the optim function from base R.
• We update ξ by numerically finding the root of the gradient of the pseudo-likelihood with respect to ξ, that is, • Finally, an orthogonalization step is necessary to comply with the orthogonality constraints on the principal component functions. The details are outlined in Section S5 next.

S5 Avoiding the Use of Orthogonal B-splines
In previous functional data work (for example, Zhou et al., 2008;Liang et al., 2020), orthogonal B-splines are proposed for estimation of the principal components. That is, for the response, for a functional basisB(·) = (b 1 (·), . . . ,b P (·)), to ensure that they propose the conditions We briefly remind the reader how one findsB(·) from a standard B-spline basis B(·), and afterwards show how this procedure can be avoided. Let where R is the P × P Cholesky matrix of the inner product matrix Ξ 0 . Then defining Finally, let Ξ 2 = B (2) (t) B (2) (t) dt andΞ 2 = B (2) (t) B (2) (·) dt be the penalty matrices for each of the bases. These penalty matrices satisfy Ξ 2 = RΞ 2 R , sincẽ using the linearity of the differentiation operator.
Enforcing this orthogonal constraint onB(·), however, removes much of the nice sparsity properties of the B-splines. We will now show that with a proper orthonormalization step, one need not enforce (15) and (16). For notation, letB i ∈ R n i ×P denote the matrix of basis evaluations at the indices in T for location i. We ignore weights for principal component computations, MCEM weights, and cluster membership for simplicity.
Lemma 1 (Penalized Least Squares). Consider the two penalized least squares problems which are the respective solutions under regular B-splines and orthogonal B-splines, respectively. Then the fitted values from the two problems are the same, that is, B(·)β 1 =B(·)β 1 .
Proof. Using penalized least squares, the two solutions are: By replacing matrices using the relationships listed above, where we use the fact that (AB) −1 = B −1 A −1 .
Thus, using either type of B-spline gives the same fitted values, since the bases are just rotated versions of each other. However, after estimation for multiple principal components, the resulting functions are not necessarily orthonormal, which motivates an orthonormalization step in both cases. Consider the univariate setting, and following Zhou et al. (2008), letΣ be the empirical covariance matrix of the scores. Also, let Θ be an estimated matrix of principal component functions for one group. Then the orthogonalization step is based on the eigendecomposition: where D is diagonal and O is orthogonal. We can see that ensuring that the resulting functions are orthonormal. Also, letΘ be the principal component functions under orthogonal B-splines.
Lemma 2. After estimating principal components under B-splines and orthogonal Bsplines, the orthonormalization step described above results in the same fitted functions for each basis.
Proof. Using orthogonal B-splines, one eigendecomposes In Figure S1, we plot the number of non-zero entries in the basis matrix for regular and gigabytes for orthogonal and regular B-splines, respectively, with 64 basis functions for temperature and salinity each. When using large amounts of additional Argo Core data when predicting oxygen, using regular B-splines is necessary for manageable data analysis.
The complete orthonormalization step then consists of orthonormalizing separately Θ g,X and Θ g,Y for each g as described above and as in Zhou et al. (2008), while also adjusting Λ g appropriately.

S6 Additional Simulation Plots and Results
We give further plots on our simulation studies, beginning with the first simulation. To further describe the setup of this simulation, in Figure S2, we plot the true principal components used for the first simulation. In Figure S3, we plot the effective sample size using importance sampling for the first simulation, as well as the adjusted Rand index Figure S2: Principal component functions in the first simulation. The -th PC of the g-th group is the (g + 2 )-th function of an orthogonal spline basis with dimension 11 with uniformly space knots. Turning now to the second simulation study, in Figure S4, we plot other principal component and linear transformation functions. In Table 3, we display the estimation results for marginal variance and range scale parameters for the principal component random fields of the predictor random field X .
For the SOCCOM data, we remove measurements with: a) salinity less than 28 practical salinity units (PSU), b) latitude greater than −25 degrees, c) pressure less than 0 or greater than 2,000, d) oxygen quality control denotes "bad" (8) e) temperature, salinity, or pressure quality control is not good (0) f) missing pressure or location data (this keeps linearly interpolated locations while under ice part of the data).
Following these reductions, we remove profiles with • any measurements of salinity less than 33 PSU or greater than 37.2 PSU • a pressure gap of 200 decibars or more between consecutive measurements • a minimum pressure greater than 100 decibars or a maximum pressure less than 1,450 decibars • fewer than 15 measurements • had a potential density estimated derivative less than −0.05 kg/L/decibar • there are fewer than five oxygen measurements with good (0) or probably good (4) quality control • an oxygen measurement greater than 400 µmol/kg Finally, for the response, we only consider data with good (0) or probably good (4) quality control.
For Core Argo data with temperature and salinity only, we make similar, but not identical, exclusions for profiles (note that the Argo and SOCCOM data have different quality flag meanings): • any measurements of salinity less than 33 PSU or greater than 37.2 PSU • a pressure gap of 200 decibars or more between consecutive measurements • a minimum pressure greater than 100 decibars or a maximum pressure less than 1,450 decibars • temperature, salinity, or pressure quality control is not good (3, 4) • missing pressure or location data (this keeps linearly interpolated locations while under ice part of the data) • the pressure adjusted error was greater than 16 decibars • the closest SOCCOM profile was greater than 500 kilometers away Observations were also removed with pressure less than 0 decibar or greater than 2,000 decibars.

S8 Details for the Initialization
Here we provide more details on the initialization procedure. First, we run k-means for interpolated values of the data to obtain an initial cluster assignment. After this initialization, we use the model introduced in Zhou et al. (2008) and extend it to a Gaussian mixture model. Let Y i and X i be defined as in the main document, and let Y 1 , . . . , Y n and X 1 , . . . , X n be the underlying functional random variables at each location. We now assume that the pairs (X i , Y i ) are n i.i.d. functional random variables. Furthermore, we have the unobserved i.i.d. latent variables Z 1 , . . . , Z n uniformly distributed on {1, . . . , G}.
We assume that the Y i admit the eigendecompositions and the random variables X i admit the eigendecompositions where the { il } are i.i.d. mean zero Gaussian random variables with varianceσ 2 Y , and the X k are i.i.d. Gaussian random vectors with independent components whose marginal variances areσ 2 X k . Let α i = (α i1 , . . . , α iQ 1 ) and β i = (β i1 , . . . , β iQ 1 ). Furthermore, let To model the correlation between X i and Y i , we use a linear regression model for the principal component scores, that is Using the similar notation as in the main document, we approximate The parameters are then estimated similarly to the procedure in the main document, with the only change being that the MCEM procedure is simpler as all spatial correlation is ignored so that obtaining the samples α for τ = 1, . . . , n is straightforward. Most M-step updates are analogous to the ones outlined in Section S4, and the M-step updates for {Σ g,α }, {Σ g,β } and {C g } G g=1 can be found in Zhou et al. (2008).
After estimating the above model by running the MCEM algorithm for a few iterations, we initialize as follows Then, we compute the residuals based on this model based on the last MC iteration. That is we computeẼ and assume that there exist functions ϕ 1 , . . . , ϕ Q 2 and independent random fields Again, we approximate B(t)Θ g,E = (ϕ 1 (t), . . . , ϕ Q 2 (t)).
We then run a single M-step update forΘ g,E and initialize Θ g,E =Θ g,E .

S9 Description of Prediction Formulas
As in other parts of the supplement, we ignore the dependence on clusters in the following developments. In Proposition 1 and Section S2, we found that, in each cluster, where X 0 and Y 0 are centered versions of the data, A =    Σ α 0 0 Σ η    , and U defined in Section S2.
We now consider new scoresα andη where we want to predict, potentially with additional predictor dataX. For simplicity, suppose thatX is available; however, we discuss later how to form predictions whenX is not available. Consider the matrix where Σ α,α = Cov(α,α), Σα = Var(α), and likewise for the other matrices involving η.
Also, note thatǍ −1 can be formed using a Vecchia approximation for the entire matrix Recall the definition of U in Section S2, and similarly define the matrices where UX ,α = D −1/2 X BXΘ X is the respective matrix for the new data. The three groups of rows correspond to X and Y , andX, respectively, while the four groups of columns correspond to α,α, η, andη.
Proposition 3. In this case, the conditional distribution of the scores can be written: where X 0 , Y 0 , andX 0 are centered versions of the data.
Proof of Proposition 3. We sketch the proof, which is similar to that of Proposition 2.
Define DX analogously to D X andB to be the basis matrix forX. Furthermore, define Each entry of 0 added to the matrix represents the scores at the new location. This gives the conditional covariance matrix Σ a , and following the steps of the proof of Proposition 2 directly gives the conditional mean as well.
This shows that prediction is quite similar to the computations in the E-step. Likewise, whenX is not available, we form U all with one group of rows deleted: and compute the same formulas with only X and Y as the available data.

S10 Product
The predicted clusters, oxygen, temperature, and salinity values, as well as their uncertainties, are given and described in the Google Drive folder https://drive.google.com/ drive/folders/1_rG5woV45UHsLpX7sWl3C-lWRaNvZysX?usp=sharing.
In Figure S5, we plot the predicted clusters for each of the 12 months.

S11 Additional Data Analysis Plots
We show additional data analysis plots of interest in a variety of figures. In Figure S6, we plot our AIC-based model selection criterion for different number of principal components.
In Figures S7 and S8, we demonstrate the seasonality dependence of the cluster memberships, showing fairly different clusters for January versus July. In Figures S9, S10, and S11, we plot predictions, predictive variances, and means for temperature and salinity, respectively, which are obtained naturally from our analysis of temperature, salinity, and oxygen.
Similarly, we plot the first principal component of each estimated in Figure S12: principal components of temperature and salinity, linear transformation functions, and principal components for E. In Figure S13, we break down the oxygen variability by the cluster. We plot the length of uncertainty estimates in Figure S14 in the leave-out experiments, showing that, as expected, predicted uncertainties for leave-float-out are higher than for leave-profile-out. In Figure S15, we give an account of the proportion of variance explained in oxygen by location, and in Figure S16 we consider the two contributions to the oxygen uncertainty in (10) and (11). In Figure S16 we can identify areas where the cluster membership is uncertain. In Figure S17, we plot predicted functions over varying latitude with fixed longitude. In Figure S18, we plot the estimated log-likelihood (as calculated in Section 3.4), as a function of the EM iteration, showing that our estimation strategy estimates parameters that increase the log-likelihood. In Figure S19, we plot predictions and variances for nitrate. Figure S6: Akaike Information Criterion values for varying number of principal components for the covariates and the response. In general, the AIC seems to be decreasing as a function of both principal components, though it increases from 12 to 13 response PCs. Figure S7: Marginal membership information for January (left, summer) and July (right, winter). This demonstrates seasonal differences in oceanographic fronts.     Figure S13: Breakdown of oxygen variability by cluster based on the model). Measurement error standard deviation refers to σ Y , standard deviation where temperature and salinity is known is σ 2 Y + B(p)Θ Y,g Var(η g )Θ Y,g B(p) , while the total standard deviation is σ 2 Y + B(p)Λ g Var(α g )Λ g B(p) + B(p)Θ Y,g Var(η g )Θ Y,g B(p) . Figure S14: The median predicted standard deviations for leave-out experiments, by pressure. Figure S15: Proportion of oxygen variance explained by the covariates and spatial correlation in prediction, January 15, 2020.  (11). The combination of this uncertainty represents the overall uncertainty in the oxygen estimates. Figure S17: Estimated oxygen (left) and oxygen standard deviation (right) for longitude -90 • for varying latitude and pressure. As one moves from south (dark) to north (light), the predicted oxygen and their uncertainties take different shapes. Figure S18: Log-likelihood f (X, Y ) for each MC iteration when considering G = 5, Q 1 = Q 2 = 12.