Water column studies focusing on the carbonate system in the GOM shelf and upper slope (up to 1000 m) region include GOMECC cruises 1–4 in 2007, 2012, 2017, and 2021 (Barbero et al., 2024, 2017; Peng and Langdon, 2007; Wanninkhof et al., 2012) ( Table 1 ). In addition to these large-scale expeditions, regional datasets also exist from individual studies focusing largely on the northern GOM shelf (Cai et al., 2020, 2011; Hu et al., 2017; Huang et al., 2015).

From April 2021 to February 2023, five cruises were conducted in the nwGOM shelf and upper slope to collect both surface and water column data across different seasons as a part of the project “Ocean Acidification at a Crossroad” (XR) ( Figure 1 ). Among the stations sampled were those along the Galveston transect, which is a transect along the 95° W longitude first taken during GOMECC-1 cruise. This transect was sampled on GOMECCs 3 and 4 and again by the XR project in April and August of 2021, December of 2022, and February of 2023 ( Table 1 ; Figure 1 ). The combined dataset before circulations QC and preparation included 960 data points from 36 stations, out of which 94 of these data points from the GOMECC-4 cruise were to be used for testing the models while the remaining 866 data points were to be split (90:10) and used for model training and calibration.

Seawater sampling for the XR cruises was done following the best practices for carbonate chemistry analysis (Dickson et al., 2007). Briefly, samples were taken from the Niskin bottles into 250-ml ground-neck borosilicate glass bottles. Preservation of samples for carbonate chemistry analyses was done by adding 100 μl of saturated HgCl₂ solution. Glass stoppers with the aid of Apiezon^® L grease and rubber bands were used to seal the samples. The GOMECC samples were analyzed on board the ship, including DO, DIC, and TA; nutrient samples were analyzed at the Atlantic Oceanographic and Meteorological Laboratory (AOML) (Barbero et al., 2017, 2024; Peng and Langdon, 2007; Wanninkhof et al., 2012). All XR carbonate chemistry (DIC, pH, and TA) samples were analyzed in our lab at Texas A&M University-Corpus Christi (Hu et al., 2018), and DO was analyzed on board the ship. Nutrient samples were analyzed at the Geochemical and Environmental Research Group at Texas A&M University.

For the GOMECC samples, DO was determined using an automated oxygen titrator with amperometric end-point detection (Culberson and Huang, 1987). DIC was determined using gas calibrations and Certified Reference Material (CRM) stability checks to ensure proper performance (with precisions of ± 1.37 µmol/kg) (Johnson et al., 1985). For the XR samples, DO was determined on a selected subset of the samples using Winkler titration (Culberson, 1991; Winkler, 1888), as these water samples were used to verify the DO sensor on the onboard the CTD (DiMarco et al., 2012). DIC was determined using infrared methodology on a DIC analyzer (Apollo SciTech Inc.) with CRM to ensure optimal instrument performance (with precisions of ±0.1%) (Chen et al., 2015). For all cruises, TA was analyzed using open-cell Gran titration (Gran, 1952), and pH was analyzed using a spectrophotometric method with purified m-cresol purple (Liu et al., 2011). Due to the COVID-19-caused interruption to CRM supply in early 2021, XR1 samples were analyzed using CRM Batch #181 and a homemade secondary reference (based on CRM Batch #181) taken from filtered (0.2 µm) and HgCl₂-dosed Aransas Ship Channel water. XR2 samples used CRM Batch #194, 197, and 198; XR3 samples used CRM Batch #202 and 204; and XR4 samples used Batch #204 only, while XR5 samples used Batch #204 and 205.

Speciation calculations for the GOMECC-1 data used TA, DIC, and nutrients as input variables as direct pH measurement using the spectrophotometric method was not available. Estimated uncertainties using the average surface water conditions for pH and Ω_Arag were ±0.01 and ±0.17, respectively (Orr et al., 2018); then DIC and pH were used for speciation calculations in all other GOMECC datasets as well as the XR ones without nutrient data, resulting in an uncertainty in Ω_Arag of ±0.18. Though some calculations were done without nutrients, uncertainty contributed to the combined standard uncertainty are negligible even in high-nutrient regions, except in high-nutrient waters with TA as a member of the input pair (Orr et al., 2018). Carbonate speciation calculation was done using the MatLab version program CO2SYS (Sharp et al., 2023; van Heuven et al., 2011). Carbonate dissociation constants from Mehrbach et al. (1973) refit by Dickson and Millero (1987), the dissociation constant of bisulfate reported in Dickson et al. (1990), the total boron concentration provided in Uppström (1974), and the aragonite solubility constant from Mucci (1983) were used in these calculations. As nutrient data were not available for the XR cruises due to technical issues encountered during the pandemic, we were not able to do a strict internal consistency analysis because using TA as one of the input variables would require nutrient information. However, the offset between calculated (using DIC/pH without nutrients, except for GOMECC-1) and measured TA values, was compared across all cruises to ensure that data quality was consistent among these trips ( Figure 2 ) (Patsavas et al., 2015).

When modeling carbonate chemistry, it is crucial to focus on a geographic area where local processes and water masses exert similar control. This ensures that the relationships between parameters remain consistent throughout the area and are accurately captured in the models. Hence, data used for our model development were limited spatially within the range of 27.1–29.0˚N and 89–95.1˚W. Removal of the upper water column of various depths is regularly done in some capacity in carbonate modeling studies in an attempt to remove large, short-term, unconstrained variabilities, such as surface water gas exchange, physical, biological, and seasonal changes (Juranek et al., 2009; Kim et al., 2010; Lima et al., 2023; McGarry et al., 2021). Depth ranges of the mixed layer to be removed vary by location and season. Some models do not require exclusion of air-surface water interactions due to predictor variables being unaffected by this exchange, and others have used removal of additional surface water as a means of removing the effect of seasonality (Juranek et al., 2009; Kim et al., 2010; Lima et al., 2023; McGarry et al., 2021; Montes et al., 2016). Since the models produced in this study could contain DO, removing air-sea interaction was necessary. However, the removal of the mixed layer to eliminate seasonality as done by Kim et al. (2010) would equate to a removal of more than 100 m to account for the seasonal changes in mixed layer depth, as winter mixed layer depth in nwGOM can be more than 100 m (Muller-Karger et al., 2015). To retain as much data as possible while removing seasonality, removal of the upper 30 m of data was arbitrarily chosen. Note this removal should more than cover the depth of the mixed layer in summer but is shallower than those depths in winter and spring. The predictor parameters (salinity, temperature, pressure, and DO) were normalized to avoid collinearity, which causes computational problems that make parameter estimates unstable (Quinn and Keough, 2002). This was done through computing Z-scores using Equation 3, where X_i is the predictor variable data point, X_m is the mean of the calibration set, and X_SD is the standard deviation of the calibration set. This can be described as centering by subtracting the mean and dividing by the standard deviation, resulting in predictor variables having a mean of 0 and a standard deviation of 1.

This approach uses normalized predictor variables to estimate nonnormalized target variables, so the relative importance of multiple parameters can be compared within one model through the magnitude of the coefficients. The normalized predictor variables were tested for collinearity by calculating the variance inflation factor (VIF) for each predictor. Combinations of terms resulting in a VIF greater than 10, indicating excessive collinearity, were not used in any model (Kutner et al., 2005). Outliers were identified using studentized residuals, which are more effective for detecting outlying Y observations than standardized residuals. Studentized residuals larger than 3 (which are classified as outliers) were further examined. Each data point was individually evaluated, considering corresponding measurements from the same sample for consistency. Data flagged as questionable across multiple parameters were carefully considered, and if multiple indicators suggested potential error, the observation was excluded from the dataset. Any data points with studentized residuals >5 were removed due to the likelihood of being erroneous or anomalous.

The modeling exercise was conducted in Python using the Scikit-learn packages (Pederegosa et al., 2011). Using a combination of the measured data (salinity or S; temperature or T; DO; and pressure, or P), these techniques produced RF and MLR models for pH and Ω_Arag. The predictor variables were selected by dredge for MLR models or variable importance for RF models (see below for details). The process was repeated until either R ² values were <0.90 or only one predictor variable remained. Final predictor variable selections were used to create final models that were evaluated.

MLR - Empirical models were developed for each target variable (i.e., pH and Ω_Arag) using least squares multiple linear regression on each combination of hydrographic variables, following the methods of Juranek et al. (2009) and Alin et al. (2012). All parameters (DO, S, T, P) in the training data were given as predictors to a full MLR model. Dredge methodology was used to evaluate every possible combination of the four parameters. From among the combinations, all models were ranked based on minimizing the Akaike information criterion (AIC), a metric for model selection that incorporates both the goodness of fit and a penalty for complexity (Burnham and Anderson, 2004). The best performing models with R ² > 0.90 and correlation coefficients with correct signs (±) were selected for further evaluation.

RF - An initial RF model was created including all four variables. Hyperparameter tuning is the process of optimizing a model’s configuration settings to improve performance by specifying tuned values for the number of trees, max tree depth, and minimum samples split. Hyperparameters for tuning were decided using scickit-learn GridSearchCV ( Supplementary Table S2 ). Selected hyperparameters were applied using an RF regressor function (Probst et al., 2019). After this step, the variable importance was analyzed, and a stepwise procedure was followed to create a new model without the variable of least importance. The new model was then tuned and evaluated, if the model's R2 > 0.90 it would be considered in the final model valuation. In addition, the process starting by analyzing feature importance to create a simpler model was continued until the resulting model did not meet previously stated criteria.

The models were used to predict their target variable while using predictor variable data from the completely “unseen data,” or data uninvolved in the training process, in this instance, the GOMECC-4 data. By using data from a cruise, which also took place within the time period that the modeled data were collected, separate from the training data, this ensures that the model performance results will accurately represent the models’ efficacy on new data. To ensure that the test dataset is representative of the training dataset and that the distributions of the variables are not problematically different, a two-step evaluation process was conducted. First, the distribution of each variable in both the training and test datasets was visually examined using graphical representations, such as histograms, density plots, or boxplots, to identify any noticeable differences. Second, the Kolmogorov–Smirnov (KS) test, a non-parametric statistical test, was applied to each variable to evaluate whether two samples are drawn from the same underlying distribution by comparing their cumulative distribution functions (Massey, 1951). These analyses collectively concluded that the training and test datasets are statistically comparable, minimizing the risk of bias in model evaluation due to distributional differences.

After confirmation of the representative quality of the test data, model evaluations were based on the performance across multiple metrics, including mean absolute error (MAE), mean square error (MSE), mean absolute percent error (MAPE), correlation coefficient R², and spatial bias (Vujovic, 2021). MAE and MSE represent the absolute-difference and the squared-difference metrics, respectively. MLR model variable coefficients were used for interpretation of relative variable effect. The RF model variable importance best describes model fit, not direct variable relationships, and therefore was not interpreted this way, that is, variable effect. Additional analysis of the residual data to identify normality and potential signals was done using QQ-plots and Shapiro–Wilks tests (Chen et al., 2019) ( Supplementary Figures S1 - S8 ).

The nwGOM area (latitudes 27.1–29.0˚N and longitudes 89–95.1˚W) was selected to isolate a relatively small geographic area where controlling factors on the carbonate system are expected to remain similar in the examined window of time and included the same stations surveyed repeatedly along the same transect (i.e., the 90° line). This study area lies almost entirely on the Texas shelf and far enough to the west so that the effects of the freshwater discharge transported westward are less pronounced than on the bulk of the Louisiana Shelf (Androulidakis et al., 2015; Morey et al., 2003). Under heavy influence from the large river’s watershed, salinity would be a key variable for models, although it is not detected as such in our models. Therefore, salinity was only included as a predictor variable in one of the Ω_Arag models and only half of the pH models.

Removal of the surface layer of the water column was done strategically to remove large, short-term temporally unconstrained variabilities, such as surface water gas exchange, physical, biological, and seasonal changes (Juranek et al., 2009; Kim et al., 2010; McGarry et al., 2021). In the nwGOM, removal of shallow water can potentially mitigate the minimal effects of freshwater outflow as the low-salinity river discharge is buoyant. Removal of surface water can be less important when air-surface water interactions do not affect predictor variables. Even still, removal can serve as a means of removing the effect of seasonality, though depth ranges of the mixed layer to be removed vary by location and season (Juranek et al., 2009; Kim et al., 2010; McGarry et al., 2021; Montes et al., 2016). Since most of the models produced in this study contain DO, removing air-sea interaction was necessary. Based on the final models chosen for Ω_Arag and pH, the relationships showed that the exclusion of the upper 30 m of the water column did improve the performance of all Ω_Arag models. When models are compared with the model containing the same spatial data but all depths, models with the upper 30 m removed, except for pH MLR models, had either greater R², lower MSE or MAPE, or a combination of the three ( Supplementary Table S2 ). The improvements in all models indicate that the removal of surface water increases the performance of the models, though it should be noted that there may be bias across seasons due to unaccounted-for surface water and mixed-layer inclusion in the datasets used for training. Seasonal bias could also occur due to a greater quantity of summer observations in the training dataset. Further data beyond the 2007–2023 time series, as well as the inclusion of more seasonal data, are needed to properly assess the potential existence of these biases.

There are many known relationships between Ω_Arag and physical and chemical parameters, which is why the variables salinity, temperature, pressure, and DO were chosen in developing Ω_Arag models. However, the degree to which each of these variables is intertwined with Ω_Arag is not so easy to identify but can be inferred by their selection and importance in the models with the highest performance.

Temperature was the most prominent predictor for Ω_Arag in all models ( Table 2 ). The strength of the relationship can be explained by Ω_Arag as a function of K_sp and [CO₃ ²⁻], and K_sp as a function of temperature, salinity, and pressure. Temperature as a dominant predictor variable aligns well with models built for Ω_Arag prediction in other areas, which also found temperature to be a primary predictor for Ω_Arag (Alin et al., 2012; Juranek et al., 2009; Kim et al., 2010; McGarry et al., 2021). Pressure was included in all MLR models and two out of three RF models ( Table 2 ). The pressure term in the models can be in part explained by Ω_Arag as a function of K_sp and [CO₃ ²⁻], and K_sp is also a function of pressure, that is, increasing solubility with depth. DO was included in both MLR models and one RF model ( Table 2 ). The inclusion of DO is likely due to the close coupling of DO and DIC in subsurface waters, which includes CO₂ uptake (photosynthesis) and release (respiration) and lacks exchange with the atmosphere in the subsurface (>30 m). Although salinity was also tested, it was only included in one MLR model based on selection criteria and was the least important predictor variable by coefficient ( Table 2 ). The lack of salinity could be in part because salinity in the modeled dataset has small variability from 33.2 to 36.7 without outliers and an interquartile range (IQR) of from 35.2 to 36.4. This means salinity has a relatively small contribution to variability, either due to the hydrography of the area or possibly in part to the limitation of the dataset’s spatial and depth coverage. By focusing on the nwGOM and removing the upper 30 m it is possible that freshwater outflow transported from the Mississippi-Atchafalaya watershed westward over the Louisiana-Texas shelf in the surface water has been partially or wholly avoided (Morey et al., 2003).

Although anthropogenic CO₂ is the dominant driver for long‐term, that is, multidecadal, change in pH in the upper ocean water column, the variations observed in coastal waters on decadal time scales are largely attributed to driving forces including production shift due to changes in nutrient input (Wang et al., 2013) and ocean circulation changes (e.g., upwelling) (Feely et al., 2008) and terrestrial influences (Cai, 2003; Duarte et al., 2013; Gomez et al., 2021). The fact that DO was included as an important predictor variable suggests that biological activities are a major driving force for the carbonate system. Similarly, DO was also found to be the strongest predictor of pH in the Northeast US (McGarry et al., 2021). Pressure was also important in all the models, as carbonate system equilibria are affected by pressure ( Table 2 ). There is a visually notable trend of overprediction of pH from 400 to 800 m; this is an area with lower pH and Ω than the rest of the water column. We attribute this bias to the non-linearity of the carbonate system. This type of bias has also been seen in other similar models in different geographic locations as well (Juranek et al., 2011, 2009).

One of the primary motivations behind the use of the MLR modeling with standardized variables is the ability to use empirical relationships among predictor variables to accurately describe the controlling processes in the study area, whereas RF modeling can only provide relative feature importance, which is not a quantification of the correlation between variables but of the relative impact on the performance of the model, or reduced impurity, that can be attributed to that predictor variable (Breiman, 2001). Due to standardization, the strength of the predictor variable correlation with the target variable can be inferred by the absolute value of the coefficients of predictor variables included in the models. However, it is important to note that statistical models reflect correlations between predictor variables and target variables but do not necessarily reveal mechanistic relationships (Quinn and Keough, 2002) ( Table 2 ).

Because of differences in hydrography, river influences, and location of the drivers of the carbonate system, statistical models for carbonate system parameters are often specific to the geographic areas where they were created for. Although many of the primary drivers of the carbonate system remain the same, the secondary variables shift, as does the degree to which they influence.

In the literature, Ω_Arag models have been created for the northwestern Atlantic (McGarry et al., 2021), northern Gulf of Alaska (Evans et al., 2013), Central Oregon Coast (Juranek et al., 2011, 2009), southern California Current system (Alin et al., 2012), Queen Charlotte Sound (Hare et al., 2025), U.S. east coast (Li et al., 2022), and the Sea of Japan (East Sea) (Kim et al., 2010). These models all have between two and three variables and an adjusted R ² ≥ 0.91 compared to between one and four variables with an adjusted R² ≥ 0.94 in this study. The primary predictor variable in all previous models except two (Evans et al., 2013; Juranek et al., 2011) is temperature, as is the case in this study. Most models also include some combination of temperature with salinity, oxygen, and interaction terms; one model also used pressure (Kim et al., 2010), and one used NO₃ ^- (Evans et al., 2013). The models created here and those created previously for Ω_Arag agreed on the use of temperature as a primary predictor variable in combination with other predictor variables.

pH models have been created for the eastern US coast (McGarry et al., 2021; Li et al., 2022), the northeast Pacific (Juranek et al., 2011), the southern California Current System (Alin et al., 2012), Queen Charlotte Sound (Hare et al., 2025), and the global ocean (Carter et al., 2021). All those models also used temperature and DO (Alin et al., 2012; Juranek et al., 2011; McGarry et al., 2021; Li et al., 2022), with one model also using an interaction term between DO and temperature (Alin et al., 2012), one using nutrients, salinity, and multiple interaction terms (McGarry et al., 2021) and, lastly, one using AOU and salinity (Hare et al., 2025). They included between two and seven variables and resulted in adjusted R ² ≥ 0.89 in comparison to our models, which have between three and four variables and adjusted R ² ≥ 0.93. Our MLR models are consistent with previous models that DO was an important predictor variable for pH.

The GOMECC-1 survey in 2007 was among the first to collect comprehensive, high-quality, measurements of the inorganic carbonate system parameters in the GOM. However, the capacity for autonomous observation of ocean carbonate chemistry is growing rapidly. By combining real-time, high-resolution data with modeling methods, it is possible to create comprehensive data coverage while simultaneously providing data validation via comparison between different data sources. With the use of MLR in combination with a single Argo profile float containing DO and temperature sensors, Juranek et al. (2011) were able to create a 14-month comprehensive time series containing Ω_Arag and pH, which were verified by comparisons to independently measured values. Evans et al. (2013) applied models to field data from glider flight and a GLOBEC mesoscale SeaSoar survey (Cowles, 2002) to create complete datasets, allowing for identification of variability of Ω_Arag. Hare et al. (2025) also had success with this in Queen Charlotte Sound, British Columbia, where regression models were applied to regional autonomous glider data. Although autonomous measurements for biogeochemical parameters in the GOM have just begun (Osborne et al., 2024), successful combination of models and autonomous data collection should greatly improve geochemical data coverage and advance carbonate chemistry studies.

In addition to long-term changes, there is strong motivation to apply empirical models to year-round hydrographic data to reconstruct the seasonal cycle of carbonate chemistry. For example, Juranek et al. (2009) created a model using data collected in spring and successfully used it to predict the seasonal changes in Ω_Arag under the assumption that the seasonal variability of the primary independent variables (temperature and DO) would be comparable to the variations encountered spatially during the initial data collection. One method to mitigate effects of the lack of seasonal variation data on the relationships between physical and chemical data is by excluding depths influenced by local meteorological conditions from the analysis (Kim et al., 2010). Most models were calibrated using mostly summer data and none with subsurface seasonal data for the GOM. Our models contain data from all seasons with slightly more data points from summer months. Seasonal application is expected to be viable in the produced models for pH as it was explained primarily by biologically driven changes, that is, the coupling between DIC and DO. Seasonal application is also expected to be viable for Ω_Arag at depths of 30–300 m, as biological changes are primarily driven by DIC rather than TA (Anglès et al., 2019; Hu et al., 2018). It is expected that changes in DIC should be captured in these models due to their relationship with DO in the subsurface (Anderson and Sarmiento, 1994; Hales et al., 2005).

Furthermore, there is also motivation to apply statistical models to forecast the carbonate system conditions as OA progresses. Previously created MLRs have been applied to generate hindcasts (up to 40 years) on carbonate chemistry (Evans et al., 2013; Kim et al., 2010). However, the temporal applicability of a model is expected to be limited. To use the model for predictions beyond the timeframe of training data, even with relatively unchanged circulation and watershed conditions, it is necessary to make corrections for shifts caused by changes in the anthropogenic CO₂ inventory in the water column, increasing error beyond model uncertainty. This is because empirical models do not account for anthropogenic addition of CO₂, which changes the ratio of DIC/DO. Due to atmospheric CO₂ input, previous modeling studies have calculated that the statistical relationships will need to be updated every 5–10 years to account for changes in target variables (Alin et al., 2012; Juranek et al., 2011, 2009). The models can be used for estimation of data within the timeframe of the training data (2007–2023), with the requirement that seasonality or waterbody changes are captured, hindcasting up to 10 years (1997–2007), and predicting up to 10 years (2023–2033) (Juranek et al., 2011). For example, Juranek et al. (2009) created a model using data collected in spring, and they successfully used it to predict the seasonal changes in Ω_Arag. The assumption of this study was that the seasonal variability of the primary independent variables (temperature and DO) would be comparable to the variation encountered spatially during the initial data collection just a year prior. However, to extend the reconstruction further than 10 years, it would be necessary to account for the progressive addition of anthropogenic CO₂ to the water column over time. For example, Kim et al. (2010) used estimates of anthropogenic CO₂ invasion to subtract anthropogenic CO₂ and apply modeling estimates of DIC in the Sea of Japan (East Sea) to create a 40-year reconstruction. In our study, the total uncertainty values are 0.12 for Ω_Arag and 0.0173 for pH. The rate of pH decline at the surface during the 2013–2022 period is 0.0025 ± 0.0005 year¹⁻, while there is no discernable trend in Ω_Arag (Hu, unpublished data). Based on this, it would take 14 years for the accumulation of anthropogenic carbon to surpass two times the uncertainty for pH (i.e., the time of emergence, Turk et al., 2019). Due to the lack of trend in Ω_Arag, these models could be applicable for a longer period in this region. However, we recommend conservatively applying caution when using beyond the 10-year window. For all models but particularly Ω_Arag, it seems likely that other factors such as changing climatology may be the first to influence model applicability rather than anthropogenic CO₂ intrusion.

Technological advancements (e.g., satellites, floats, and gliders) have made autonomous measurements for selected variables possible. Even though satellites provide large volumes of invaluable data, their scope is limited not only to the surface layer but also to measurements of temperature and ocean color, which can be used for estimations of productivity, dissolved organic matter, and sea surface CO₂ partial pressure (e.g., Chen et al., 2019). Some remedies for the lack of data beyond the surface have come in the form of novel technologies, including Argo floats and gliders (Osborne et al., 2024). The Argo array supports more than 3,000 floats across the global ocean, supplying profiles of T and S, of which about 200 are also equipped with a DO sensor. Biogeochemical-Argo (BGC-Argo) floats can take measurements from the surface to 2,000 dbar and provide high-resolution autonomous measurements for some combinations of chlorophyll, particle backscatter, DO, nitrate, pH, or irradiance (Roemmich et al., 2019). Nevertheless, this equipment is limited in coverage ability; BGC-Argos have yet to reach the comprehensive coverage of traditional Argos, and they usually do not operate on continental shelves. Gliders share the similar biogeochemical sensors used by the BGC-Argo floats with a similar mode of operation but perform sawtooth trajectories from the surface to the bottom, or to 200–1,000 m depth, and their deployments are limited largely by biofouling, battery life, and payload space, hindering their ability to carry many sensors (Saba et al., 2019; Testor et al., 2019). Additionally, all data from floats and gliders must undergo rigorous quality control and correction of potential errors in sensors, such as drift and be compared to lab-measured samples when possible (Roemmich et al., 2019).

Given the performance of the models created here, it would be reasonable to expect that in combination with datasets from data collection via these platforms, one or both of pH and Ω_Arag could be estimated. Traditional Argos, as well as most gliders, have temperature and pressure (or depth) measurements, meaning Ω_Arag could be estimated using either of the two RF models with R ² values of 0.94 and 0.97. With the addition of DO, which is included in a subset of the Argos and in the BGC-Argos, the Ω_Arag model with an R ² of 0.98 could be applied. Additionally, the pH model with R ² of 0.94 could be used; even though pH sensors are included in many BGC-Argos, the modeled results could provide robust data quality check capability.

The models developed in this study serve as valuable tools for reconstructing pH and Ω_Arag data in the absence of direct chemical observations, leveraging available hydrographic information. The models can also be used for hindcasting over a 10-year period (1997–2007) and for forecasting over the next ~10 years (2023–2033) in the nwGOM, if seasonality and watermass changes are adequately captured. Through modeling, our results showed that RF models do not have any distinct advantage over the MLR models which are commonly used. It is crucial to examine potential shifts in circulation, water mass composition, and anthropogenic CO₂ accumulation to refine and update the models over time. Despite their robustness, several potential issues remain across all models, suggesting room for improvement, such as difficulties estimating target variables due to significant variation in the shallow and mixed layers, a need for additional data points, or the application of more complex models. Future work to enhance these models could involve incorporating additional data, as well as exploring the application of neural network AI models and locally interpolated regressions, such as those employed by Carter et al. (2021).