{'fraction_non_alphanumeric': 0.03806444884763792, 'fraction_numerical': 0.015469113915250113, 'mean_word_length': 4.682006715386135, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 4, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
3D Bivariate Spatial Modelling of Argo Ocean Temperature and Salinity Profiles

Variables contained within the global oceans can detect and reveal the effects of the warming climate as the oceans absorb huge amounts of solar energy. Hence, information regarding the joint spatial distribution of ocean variables is critical for climate monitoring. In this paper, we investigate the spatial correlation structure between ocean temperature and salinity using data harvested from the Argo program and construct a model to capture their bivariate spatial dependence from the surface to the ocean's interior. We develop a flexible class of multivariate nonstationary covariance models defined in 3-dimensional (3D) space (longitude x latitude x depth) that allows for the variances and correlation to change along the vertical pressure dimension. These models are able to describe the joint spatial distribution of the two variables while incorporating the underlying vertical structure of the ocean. We demonstrate that proposed cross-covariance models describe the complex vertical cross-covariance structure well, while existing cross-covariance models including bivariate Mat\'{e}rn models poorly fit empirical cross-covariance structure. Furthermore, the results show that using one more variable significantly enhances the prediction of the other variable and that the estimated spatial dependence structures are consistent with the ocean stratification.


Introduction
The international scientific program named Array for Real-time Geostrophic Oceanography (ARGO) was launched in the early 2000s as a response to the call of global observation networks to monitor the climate system (Argo, 2000;Johnson et al., 2022). Since then, the program has launched a global network of 4000 free-drifting Argo profiling floats that measure ocean variables in the upper 2000 meters of the world's oceans. Each float performs a ten day duration "park-and-profile" mission. From the surface, the float descends to the drift depth at 1000 meters where it will park for 9 days. At the 10th day, the float descends to 2000 meters and collects temperature (in degree Celsius or • C), salinity (in practical salinity unit or PSU), and pressure (in decibars or dbar) measurements as it ascends to the surface. Once at the surface, the data collected are transmitted via satellite. The full "park-and-profile" mission is illustrated in Figure 1. The data harvested by the floats are made available within 24 hours of its collection as Argo data products (https://argo.ucsd.edu). By making the data publicly accessible, the Argo program offers the research community the opportunity to analyze ocean processes and challenges modelers to improve existing data-driven scientific methods.
The core Argo data products are temperature and salinity measurements, the two most important variables through which other oceanographic variables, such as freezing point, electrical conductivity, viscosity (Pawlowicz, 2013), ocean heat content and potential density (Yarger et al., 2022), and tropical cyclone-induced ocean thermal response (Hu et al., 2020), can be derived. The temperature and salinity variables have been key pieces in understanding the physical properties and dynamics of the ocean. Their distributions were shown to drive ocean circulation Gangopadhyay, 2022), affect climate processes (Olson et al., 2022), and change biogeochemistry (Galán et al., 2021;Ding et al., 2022). They were also used to understand underwater sound propagation for hydroacoustics research (Jana et al., 2022), design offshore wind farms (Escobar et al., 2016), and identify stressors to organisms in seawater and freshwater (Walker et al., 2020). Figure 1: Illustration of the standard Argo "park-and-profile" mission. Source: (Wong et al., 2020)

Spatial Interpolation of Argo Data
Spatially continuous maps of temperature and salinity measurements are essential to multidisciplinary scientific research. However, despite the planetary scale and subsurface reach of the Argo profiling network, not all locations are sampled. Since other ocean variables depend heavily on temperature and salinity, efforts are focused on obtaining the best interpolated maps of these two variables. Several institutions have produced high-resolution gridded global temperature and salinity datasets from the sparse measurements recorded by the Argo profiling floats using various interpolation techniques (Liu et al., 2020). The list includes the EN4 dataset from the U.K.
Met Office (Good et al., 2013) and the grid point value of the monthly objective analysis using Argo data (MOAA) of the Japan Agency for Marine-Earth Science and Technology (JAMSTEC) (Hosoda et al., 2008), both of which use various optimal interpolation methods and covariance functions. Another example is the Barnes objective analysis Argo (BOA-Argo) gridded dataset that employs objective interpolation based on Barnes successive correction method (Li et al., 2017).
One of the most highly used Argo data product is the Roemmich-Gilson Argo climatology from the Scripps Institution of Oceanography (Roemmich and Gilson, 2009). They used weighted local regression in its high-resolution global mapping. Suppose s = (L, l, p) is the spatial location vector and t is time in yeardays. For s 0 = (L 0 , l 0 , p 0 ) , a reference spatial location vector, and t 0 , the reference time index, they fit the mean function: (1) to measurements from the sampled location (L, l, p), such that (L, l) is one of the 100 nearest neighbor from (L 0 , l 0 ), for t within 12 calendar months from t 0 , and for p which is one pressure level above and below p 0 . Here β 0 and β k , γ k , and δ k , k = 1, . . . , 6, are scalar coefficients.
Previous studies begin with discretizing the vertical domain of the Argo data into distinct pressure levels and performing the high resolution spatial interpolation of each layer separately.
Furthermore, prior to fitting their proposed means and covariance functions, given the irregular sampling of Argo profiles along pressure, a linear interpolation is done to obtain measurements at the fixed pressure levels. This extra step can significantly introduce errors to the model. To address this problem, Yarger et al. (2022) proposed treating the Argo profile data as a functional data that changes with pressure and introduced a functional model for the mean of the form: Here y denotes the year the measurement was obtained, 1 (·) is an indicator function, and β 0,y (p) and β k (p), k = 1, . . . , 6, are functions that direct how the coefficients change along pressure.
Similar to the mean function in (1), the functional mean (4) is fitted to the measurements from the nearest neighbors. This functional data approach reports lower prediction errors than the pressure-by-pressure approach. The residuals from Yarger et al. (2022) were then modelled using functional principal components.
The aforementioned studies and other results in the Argo literature typically model temperature separately from salinity rather than leverage joint information that may be contained by the two variables. In this work, we extend the spatial modelling of temperature and salinity residuals to include their joint spatial dependence structures. Furthermore, we propose bivariate spatial models in 3D space that can incorporate special oceanographic structures. The ocean, as a fluid, has a density which is a continuous function of depth or pressure. This density dictates the layering of the ocean waters or ocean stratification which is a phenomenon that occurs when different types of water meet and mix, e.g. low-temperature and low-salinity subarctic water, high-temperature and high-salinity subtropical water, and other mixtures (Sambe and Suga, 2022). Because temperature and salinity exist within this fluid, the ocean stratification can affect their distributions. Bivariate spatial models that account for the ocean structure in the vertical dimension certainly will provide a better description of the processes involved.
The rest of the paper is structured as follows: Section 2 describes the data used for the analysis and provides various empirical analyses. Section 3 introduces the models for the 3D bivariate spatial modelling. Section 4 presents the results of model fitting and prediction. Section 5 concludes.

Background
We aim to uncover the spatial relationship lingering among the temperature and salinity residuals obtained by Yarger et al. (2022) so that a clear understanding regarding the dynamics between temperature and salinity can be made. Accomplishing such task can enable one to visualize the vertical structure of the ocean that is not directly observable and to exploit their spatial relationship for prediction at locations with no measurements. More complex patterns can be revealed as we dive below the surface. SupposeD is the area  surrounding a reference location s. Using the following formulas: we compute for the empirical variances and colocated correlation between temperature and salinity, respectively, at every location along the 175 • W and 30 • W longitude lines at a 1 • × 1 • × average, up to 100 meters down and is directly heated by the sun. This layer is where the winds, waves, and currents occur such that their interaction allows for the homogeneous dispersion of heat throughout its depth. The turbulent mixing happening in this region produces uniform vertical structure of temperature and salinity measurements, hence, the variances at this zone are generally not as high as those seen in the next zone (Helber et al., 2012;Maes and O'Kane, 2014). The pycnocline is the region where density changes drastically with depth. This layer is an invisible blanket that separates the warmer and less dense surface layer above and colder and denser deep layer below. This region coincides with the thermocline, the region of sudden decrease in temperature, and with the halocline, the region of sudden increase in salinity. A salient feature of this zone is that it is where the highest variances appear. This is because most eddies are trapped in this particular region which introduce chaotic structures (Dong et al., 2014;Chen and Wang, 2016;Chen et al., 2018). Past the pycnocline is the intermediate zone or the deep layer, a region most isolated from the atmosphere. In this layer, the temperature and salinity values are fairly uniform which explains the sustained small variances from the boundary of the pycnocline.
Figure 4 provides valuable insights on how the two processes interact. Generally, temperature decreases with pressure while salinity increases with pressure. However, this relationship does not directly translate to the observed correlations. From the plot, it can be seen that within a thin region at the surface, temperature and salinity are generally independent, with correlation more or less equal to zero. As they approach the pycnocline, they become highly positively correlated with peaks recorded at the depths of 200 to 500 meters. Beyond the pycnocline, along the 175 • W longitude, the correlation drops again as pressure increases until temperature and salinity become highly negatively correlated. This behavior of peak correlation at the pycnocline was also observed in Chen and Geng (2019). According to the authors, this pattern is not surprising as the climatological patterns of temperature and salinity are homogeneous around this area. As a result, this layer is also called the "temperature and salinity mirror layer". A similar pattern is not followed by the locations along the 30 • W longitude where high positive correlations can persist from the surface to 2000 meters down. This sustained intense correlation can be explained by the dynamic process called meridional overturning circulation which is strongly experienced in the Atlantic Ocean crossed by the 30 • W longitude (Chen and Geng, 2019). A process called density compensation allows for this temperature-salinity relationship to hold in this area.

Methodology
In the previous section, it was shown that the variances and colocated correlation of the residuals change with depth and the vertical structure at one location may not be representative of the other locations. For such phenomenon, spatial nonstationary cross-covariance functions defined in 3D space ought to be used to model the second-order structure of the residuals. We now build 3D spatial models that can accommodate the aforementioned spatially varying features, namely, the variances and colocated correlations, that are present in the Argo residuals.

Multivariate Nonstatonary Spatial Models in 3D domain
We treat temperature and salinity residuals as variables of two different 3D (latitude × longitude × depth (i.e., pressure)) processes with zero mean that may or may not be correlated and encode the information regarding their spatial dependence structure in a spatial nonstationary cross- that Z 1 = T and Z 2 = S. Here Z i (s 1 ) is the value of variable i at location s 1 and Z j (s 2 ) is the value of variable j at location s 2 . Only functions that verify the well-known requirement of positive definiteness can be considered valid models for C ij .
The literature is replete with cross-covariance function models. The parsimonious Matérn cross covariance function proposed by Gneiting et al. (2010) is one of the most popular models through which more sophisticated models are built upon. It has proposed the form: is the modified Bessel function of the second kind of order ν, Γ(·) is the gamma function, and q is the number of variables. Here the parameter list includes, for i = j, the marginal variance (σ 2 i > 0), and smoothness (ν ii > 0) parameters, and for i = j, the colocated correlation (ρ ij ) and cross smoothness (ν ij > 0) parameters, such that ν ij = 1 2 (ν ii + ν jj ) and for any σ 2 ii , ν ii > 0, d ≥ 1. Here (β ij ) q i,j=1 is a symmetric and positive definite correlation matrix.
For 3D processes on a globe such that s = (L, l, p) represents the latitude, longitude, and pressure, the following function can be used: where a h and a v are the scale parameters in the horizontal and vertical directions, respectively, and ch(L 1 , L 2 , l 1 − l 2 ) is the chordal distance given by the formula: Here R is the radius of the sphere. For global processes such as the ocean temperature and salinity, the sphere is the Earth with R =6,371 km.
The Matérn cross-covariance function in (5) with distance function in (7) is nonstationary with respect to the latitude but stationary with respect to longitude and pressure. Moreover, the model is stationary with respect to its parameters such that the variances, σ 2 i , colocated correlations, ρ ij , and the scale parameters, a h and a v , are kept constant everywhere in the sphere, for i, j = 1, . . . , q. Kleiber and Nychka (2012) proposed a nonstationary extension to the parsimonious Matérn cross-covariance function (5) in 2D (latitude × longitude) such that the variances and scales can vary in space, i.e., where D ij (s) = 1 2 {D i (s 1 ) + D j (s 2 )} is a d × d positive definite kernel matrix which controls the spatially varying anisotropy, for i, j = 1, . . . , q. Despite these innovations, a valid extension of (5) with nonstationary colocated correlation parameter, ρ ij (s), remains elusive. Turns out, the task of generalizing ρ ij to change with location is not trivial. One way to introduce nonstationarity to ρ ij is to replace β ij in (6) with a kernel smoothed nonparametric estimator, β ij (s 1 , s 2 ), as suggested by Kleiber and Nychka (2012), and in the same principle, σ 2 i in (5) can be substituted with its kernel smoothed nonparametric estimator using the following formulas: where K λ h ,λv (s,s) = exp −{ch 2 (L 1 , L 2 , l 1 − l 2 )/λ h + (p 1 − p 2 ) 2 /λ v } is a nonnegative kernel function with horizontal and vertical bandwidths, λ h and λ v , respectively, for i, j = 1, . . . , q.
This approach, however, does not guarantee that the resulting cross-covariance matrix is positive definite. This issue can be remedied by finding the closest positive definite cross-covariance matrix which is done by setting the negative eigenvalues to a very small positive number. Kleiber and Genton (2013) provided parametric forms of β ij (s 1 , s 2 ) catering to a specific phenomenon involving minimum and maximum temperature residuals over Colorado which are positively correlated at low elevations and almost independent at high elevations. Their approach hinges on the inclusion of a nugget effect. Moreover, their method is not quite appropriate for usage on the bivariate Argo dataset because the empirical colocated correlations shown in Figure 4 are not monotonic.
A class of spatial nonstationary cross-covariance functions in 2D using differential operators with respect to latitude and longitude was developed in Jun (2011). The univariate version (in space-time setting) was first developed in Jun and Stein (2007). In their work, the multivariate spatial process can be written in the form: where X k (L, l), k = 0 . . . , K, are independent zero-mean univariate spatial processes, and a i,k (L), A restriction inherent in the model above is that the variables, Z i , i = 1, . . . , q, rely on the same X k (L, l), for any k. The model can be made more flexible by replacing the univariate spatial processes, X k (L, l), with multivariate spatial processes, X k (L, l) = {X 1,k (L, l), . . . , X q,k (L, l)} ∈ R q such that Z i do not necessarily depend on the same underlying univariate process but rather on correlated univariate processes. In the following proposition, we propose an extension of (11) with this new innovation and an additional differential operator with respect to pressure: Proposition 1 Suppose X k (L, l, p), k = 0 . . . , K, are independent zero-mean multivariate isotropic spatial processes, such that each X k (L, l, p) has q elements with cross-covariance function C 0 ij and (L, l, p) denotes the spatial location in latitude, longitude, and pressure. Define a multivariate process as follows: Furthermore, suppose that the cross-covariance function C 0 ij is chosen to be the parsimonious Matérn cross-covariance function in (5). Then the multivariate process Z = (Z 1 , Z 2 , . . . , Z q ) has a matrix-valued nonstationary cross-covariance function where The proof of the proposition is rather straightforward. Refer to Appendix in Jun and Stein (2007) for the exact expressions of h r and h rs , for r, s = 1, . . . , 4.
The new formulation has two main features not present in (11). First, for any k, each X i,k (L, l, p) can have different smoothness, for i = 1, . . . , q. Second, the pair, X i,k (L, l, p) and X j,k (L, l, p), for any k, has a colocated correlation ρ ij , that may not be one for i = j. Note that the model in (11) results in ρ ij = 1 for any i, j, k.
Based on the form of the cross-covariance function in (13), it can be seen that nonstationarity can be introduced in the marginal variances, C ii (L, L, 0, p, p), and cross-covariance, C ij (L, L, 0, p, p), along the pressure dimension through the functions c i,k (L, p) and c j,k (L, p).
While the parameters a i,k (L) and b i,k (L) introduce asymmetry in the latitude and longitude dimensions, respectively, for the field Z i (L, l, p), and ρ ij is the colocated correlation between  Figure 4, respectively, when ρ 12 = 1. Figure 5 shows the corresponding colocated correlation curves. From the plots, it can be seen that the values of a i,1 (L), b i,1 (L), c i,1 (L, p), and c j,1 (L, p) influence the shape of the colocated correlation curve. The parameter ρ 12 , on the other hand, acts like the amplitude of the curve such that it controls the horizontal scaling of the curve.

Models
We pursue the modelling of the Argo residuals profile by fitting the following seven models with different features and comparing their performances at the six reference locations highlighted in Section 2. Here, the "independent" models allow different marginal structures but the correlation between the two processes (temperature and salinity) is assumed to be zero.
where M is the number of basis functions, {B i,m (p)} M i m=1 are the B-spline bases constructed using M i knots, and {β i,m } are the corresponding weights of the B-spline bases for variable i.
The number of knots are allowed to vary between variables and among the reference locations depending on the observed dependence structure of the data. e) Models I2 and B2 may produce a non-positive definite cross-covariance matrix. In such cases, we find the closest positive definite cross-covariance matrix by setting the negative eigenvalues of the cross-covariance matrix to 10 −17 .

Estimation Procedure
The nonstationary models I3, B3, and B4 are more versatile but they are more heavily parameterized than the stationary ones. Through c i,1 (p), which is parameterized using B-splines, I3, B3, and B4 can capture nonstationary vertical dependence. However, the number and location of the knots of the B-splines will dictate the goodness-of-fit of these models. Using thinly spaced knots is ideal as more detailed nonstationary dependence structure can be captured at small depth intervals. Doing so, however, will increase the dimension of the parameter vector which might cause the estimation routine to run into convergence issues. We address these problems by (I) choosing the optimal set of knots through cross-validation, (II) retrieving reliable initial parameters for the independent model, I3, and (III) using the estimated parameters of I3 as initial parameters for the bivariate models, B3 and B4. For the knot locations, we begin by placing them at the boundaries of the ocean stratification as seen in Figure 2, i.e., at 0, 100, 1000, and 2000 meters. With this preliminary set of knots, we perform the following two-step estimation procedure for I3.
Step 1. Weighted Least Squares Estimation (WLS): Jointly estimate the parameters by minimizing the objective function: where θ is a vector of parameters andĈ ij (s a , s b ) is the kernel smoothed empirical crosscovariance computed from the data as follows: Step 2. Maximum Likelihood Estimation (MLE): Jointly estimate the parameters, with the WLS estimates from Step 1 as initial values, by maximizing the Gaussian log-likelihood func-tion: where Z = {Z(s 1 ) , . . . , Z(s n ) } is the vector of all measurements, Z(s) = {T (s), S(s)}, and Σ(θ) is the nq × nq cross-covariance matrix for Z with determinant |Σ(θ)|.
The above model fitting procedure can be repeated with a different set of knots such that new knots are added at the midpoints of the previous set of knots until the value of l(θ) for I3 does not significantly improve anymore. Once the MLE estimates of I3 are obtained, we perform the MLE of B3 and B4 with the MLE estimates of I3 as initial values. Note that B3 has the same set of parameters as I3, however, B4 has one additional parameter, ρ 12 , initialized at zero. The aforementioned fitting procedure is aligned with the overall motivation of this paper, that is, to show whether bivariate modelling can improve upon modelling the variables independently.
Despite employing a local stationary framework, the size of the Argo residuals dataset can render the operations involved for MLE computationally infeasible. Inverting the covariance matrix requires O(n 2 q 2 ) memory complexity and O(n 3 q 3 ) computation complexity. On average, a 900 km radius around a reference location contains 2,000 Argo profiles. Furthermore, an Argo profile, on average, has 100 vertical measurements. This translates to a cross-covariance matrix of size 400,000 × 400,000. A high performance implementation of the models must be sought when using all those data points. We implement the proposed model on ExaGeoStat (Abdulah et al., 2018), a high performance computing software, to accelerate the operations on the Gaussian log-likelihood function.

Results
In this section, we examine the performance of the different candidate spatial covariance functions in modelling the residuals around the six reference locations. Each reference location exhibits a particular vertical dependence structure which some models may be able to represent better than the others. We measure and compare the performances in terms of how well the model   Table 2 reports the aforementioned values.
From the log-likelihood values, the bivariate models have better performances than their independent counterparts. The most flexible model that best explain the 3D bivariate dependence structure between temperature and salinity residuals appears to be B4. The log-likelihood value of B4 at all six reference locations is the highest among the seven models. Furthermore, as I3 is  an independent nested version of B4, the performance of B4 is always better than I3. However, this is not what is always seen between the I3 and B3 modelling performances. In RL1, RL2, RL4, and RL5, the independent model, I3, has higher log-likelihood values than the bivariate model, B3. This is because I3 is not nested under B3. Due to its formulation, the estimation of B3 effectively starts with a fixed parameter ρ 12 = 1, while B4 starts with ρ 12 = 0. The restrictive form of B3 does not allow for the cross-covariance structure to deviate away from the marginal covariance structures. This is where the utility of Proposition 1 appears. The additional parameter ρ 12 in B4 allows the other parameters to adjust in order to capture the true underlying cross-covariance structure. In Table 3, we summarized the estimated parameter values associated with the differential operators in the horizontal dimension. For brevity, we omit the presentation of the parameters associated with c i (p). That the parameters of B4 deviate from its simpler counterpart, B3, suggests that B3 is way too restrictive as a bivariate model. Figure 6 displays the fitted variances and colocated correlation values by the different models.  To illustrate the benefits of having a secondary variable to aid the prediction of a primary variable, we perform prediction on one variable at a time. When predicting the temperature (salinity) residuals at the Argo float which was left out for prediction, we use the temperature and salinity residuals included in modelling fitting plus the salinity (temperature) residuals at the same Argo float which was left out for prediction. While the interpretation of the results for the modelling as expressed by the log-likelihood and AIC criteria is straightforward, the prediction metric MSE does not point to one obvious winner. Except at RL5, the bivariate models report lower MSEs compared to the independent models. This demonstrates the utility of having another correlated variable present to aid prediction. To improve the prediction at RL5, we may need to prescribe a more sophisticated form for c i (p) and the cross-covariance parameters may need to be set independently from the marginals. The prediction performances at RL5 reveal the importance of correctly modelling the colocated correlation. Because we assume that we have the measurements of the secondary variable at the same location where we aim to predict the main variable, the models generally give more weight on the colocated secondary variable compared to the values of the other neighboring Argo floats. Thus, underestimating or overestimating this colocated correlation may lead to erroneous predictions.

Concluding Remarks
Argo floats play an important role in sampling the temperature and salinity measurements of the world's oceans. Despite the global reach of the Argo network of floats, not all locations are sampled. We contributed to the efforts of producing high quality temperature and salinity interpolations at unsampled locations through a bivariate modelling of the Argo residuals. Motivated by the interesting physical dynamics between temperature and salinity, we took the modelling one step further by incorporating the vertical properties of the ocean. Consistent with previous studies, we illustrated that the correlation between temperature and salinity changes with depth and successfully modelled this vertical dependence structure using a nonstationary cross-covariance function. Our proposed model is flexible enough to accommodate different vertical oceanographic regions. This study provides a framework to improve mapping of environmental variables in 3D. Further work may be pursued such as constructing new models that allow the cross-covariance parameters to be set independently from the marginals and the range parameters to change with depth. Furthermore, a 4D modelling may also be explored with time as an additional dimension.