Introduction

Ocean acidification is an unavoidable consequence of increasing atmospheric CO21,2,3,4,5. Understanding how fast the oceans acidify is critical to examining and predicting the impact of ocean acidification on the aquatic ecosystem. Invaluable decadal ocean time series measurements of carbonate chemistry parameters can now be used for this purpose. Observations show that surface ocean pH under in situ temperature (pHinsitu) has strong seasonal variations, with a high pHinsitu in winter and a low pHinsitu in summer. Long-term decreasing trends are superimposed on the regular seasonal variations6,7,8.

To investigate how fast the oceans acidify, the widely used first-order simple linear regression (SLR) method has been used to model the long-term temporal changes in pHinsitu (βpH)2,6,7,8,9,10,11,12,13. The SLR model in general form is as follows,

where pH° is the intercept, βpH is the slope and t is time in the pHinsitu vs. t plot.

The βpH, the so-called acidification rate, refers to an average rate of temporal change in pHinsitu. Conventional wisdom has it that the longer the pHinsitu time series (e.g., over a decade) the better the estimation of the βpH. Bermuda Atlantic Time Series Study (BATS) and the Hawaii Ocean Time Series (HOT) (see Fig. 1 for the station locations) are two examples of long-term time series with even sampling distribution and their βpH values fairly reflected the acidification rates due to the increased atmospheric CO2 concentrations6,7. However, in a time series with short observation years (e.g., less than a decade), the small long-term temporal change in pHinsitu due to the increase in atmospheric CO2 could be easily masked by the strong seasonal variations, especially in the case of uneven sampling distribution. The pHinsitu and seawater temperature (T) time series at the European Station for Time Series in the Ocean Canary Islands (ESTOC) between 1995–2009 are shown as examples (Fig. 2). The raw data (blue cross), when fitted by the SLR, result in a negative βpH of −1.84 ± 0.39 × 10−3 pH unit yr−1 and positive SLR-determined rate of temporal change in T (βT) of 0.023 ± 0.039 °C yr−1 (Fig. 2). The red circles represent an extreme example of uneven sampling distribution, with the sampling time gradually shifting from summer to winter. As a result, these data now would have a positive βpH of 1.61 ± 0.37 × 10−3 pH unit yr−1 and a negative βT of −0.273 ± 0.036 °C yr−1. Moreover, the slopes of the regression lines would increase with a shift in sampling time from summer to winter within a shorter observation period. This shows that βpH is highly sensitive to the temporal distribution and the duration, of time series data.

Figure 1
figure 1

Station locations of the time series, including Bermuda Atlantic Time Series Study (BATS), Carbon Retention in a Colored Ocean Project (CARIACO), European Station for Time Series in the Ocean Canary Islands (ESTOC), Hawaii Ocean Time Series (HOT), South East Asia Time Series Study (SEATS) and 137°E repeated hydrographic line (137°E).

The map was generated using Ocean Data View 4.6.2. Schlitzer, R., Ocean Data View, odv.awi.de, 2015.

Full size image
Figure 2
figure 2

Time series of pHinsitu and T at ESTOC.

The blue crosses and lines show the raw data and the regression lines. The red circles and lines show a case of selected data and the regression lines. Data taken from Gonzalez-Davila and Santana-Casiano19.

Full size image

It has been suggested that the effect of seasonal variations on βpH could largely be reduced when pH anomaly data (observed data minus seasonal average) are used2. However, this method may not be applicable to time series with short observation periods and uneven sampling distribution, which do not have sufficient data to generate effective seasonal means, resulting in errors. Furthermore, the reported low values of the coefficient of determination (R2), ranging from 0.09 to 0.55, indicate that these time series data were still poorly fitted by SLR lines2. Consequently, estimating βpH using only t introduces a deviation due to the influence of strong temporal variations. As a result, it remains unclear whether the observed rates at which surface seawater pHinsitu changes differ are due to the changing carbonate chemistry of seawater at different time series stations or to deviation using SLR.

In this study, we propose a simple statistical method for reducing deviations of βpH under the influence of strong temporal variations. We show that when pH is expressed as a linear function of t and the other variables, βpH can be expressed as a linear combination of a constant (an estimated rate of temporal change in pH) and SLR-determined rates of temporal changes in the selected variables (deviation largely due to various sampling distributions under the strong seasonal pH variability), despite complications due to different observation durations and temporal sampling distributions. Using this method, we demonstrate that βpH differs from the observations from five time series stations (the BATS, Carbon Retention in a Colored Ocean Project (CARIACO), ESTOC, HOT and the South East Asia Time Series Study (SEATS, with bottom depth >3700m)) mainly due to deviations under various sampling distributions rather than changing seawater chemistry. Average summer and winter βpH values from 31 stations along the 137°E repeated hydrographic line (137°E) likely contain the same deviation.

There are few reported acidification rates for the world’s oceans, as there are only a handful of stations with decadal time series of carbonate parameters. Time series studies with high temporal variations and few observation years can benefit from our approach to find the fair estimation in the rate of temporal change in pHinsitu.

Results

Decomposition of β pH and its physical meaning

Statistically, for any pHinsitu time series, βpH is defined as follows14,

where pHi and ti are pH and t, respectively at t = ti and is the average t.

As t can only be used to model the linear temporal change in pHinsitu, it cannot be used to model the strong seasonal variations. Suppose that other variables (e.g. T, dissolved oxygen, nutrients), x1, …, xn, are added to model the temporal variations in pHinsitu, which can be expressed as follows,

where m1 ~ mk are constants, t° is a reference point of t and are x2 ~ xk at t = t°, respectively.

Substituting Eq (3) into Eq (2), Eq (2) can be rewritten as follows,

As and , the equation above can be simplified as follows,

where are SLR-determined rates of temporal changes in x2 ~ xk, respectively (see the definition shown in Eq (2)).

Equation (4) shows that in fact, βpH can be expressed as a linear combination of a constant (m1) and rates of temporal changes in the other regression variables (). Thus, although differences in the number of observation years and temporal distribution of sampling yield different βpH and , values for change with βpH, whereas m1 remains unchanged. Statistically, if x2 ~ xk in the above equation have no long-term trends of changes, are expected to be zero. Non-zero refers to the consequence of various sampling distributions under seasonal variations in x2 ~ xk. If the long-term trends of changes in x2 ~ xk have no collinearities with that of pHinsitu, any long-term trends of changes in x2 ~ xk would not affect βpH.

Notably, m1 is defined as the amount of pHinsitu change as t changes. Therefore, if long-term changes in do not vary with that of t or pHinsitu, m1 refers to the estimated rate of temporal change in pHinsitu due to different physical and biological forcings (e.g. changes in anthropogenic CO2 concentration, productivity, microbial respiration rate, etc.), whereas refers to the deviation due to the temporal variation. Notably, the deviation due to uneven sampling distribution under strong temporal variations can be quantified using rates of temporal changes in other variables.

Verifications and criteria of additional regression variables

To verify Eq (4), additional variables need to be added to the regression model shown in Eq (3). These variables must be strong in seasonal variations and weak in long-term changes, so that they can be used to model the seasonal variations in the pHinsitu. At the same time they should have weak collinearities with t. Seawater T may serve this purpose, as pHinsitu is related to T and observations show that the T time series routinely mirrors the pHinsitu time series6,7. Furthermore, modelled results from recent studies show that the surface ocean pHinsitu has no observable change under a long-term change in T15,16.

To verify the appropriateness of using t and T in modeling pHinsitu, Eq (3) is simplified as follows,

where To is T when t = to.

Using the initial year of t° = 1988 as a reference, the multiple linear regression (MLR) model of Eq (5) is then applied to the surface ocean pHinsitu from five time series stations. The results are shown in Fig. 3 and Table 1.

Table 1 MLR coefficients of Eq (5) of seawater pHinsitu from various time series studies. The pH° and T° are determined using the SLR method when t° = 1988 (see Supplementary Table S2).
Full size table
Figure 3
figure 3

The surface seawater pHinsitu at various time series stations.

The black and red lines show, respectively, the observed values and the values modeled by Eq (5).

Full size image

Linearity between pHinsitu, t and T.

Figure 3 shows the observed pHinsitu time series and the modeled results using MLR at various time series stations. The regression results match well with the annual, interannual and long-term variations in observations. The m1 and m2 of −1.6 ± 0.1 × 10−3 pH unit yr−1 and −8.7 ± 1.4 × 10−3 pH unit °C−1, respectively, are similar among the time series stations. The R2 values are increased from 0.09 ~ 0.55 to 0.51 ~ 0.94, with all sites averaging 0.75 ± 0.17. The differences between the measured and the MLR results are at most approximately twice the pH calculation error of ±0.0062, based on the carbonate system using total alkalinity (TA) and dissolved inorganic carbon17, suggesting the validation of Eq (5). Consequently, Eq (4) can now be simplified as follows,

It is important to point out that the linearity between pHinsitu and T is a lump sum of various physical and biological forcings. Based on this linearity, the change in βpH due to various sampling distributions under the strong seasonal variability is reflected in m2βT. The fact is that, changes in pHinsitu and T owing to other factors (e.g. global warming, enhanced upwelling, etc.) may not equal to m2, resulting in residuals. The standard errors of Eq (5) shown in Table 1 possibly reflect the errors from minor disturbances and from the measurements. Since Eq (6) is generated from Eq (5), Eq (6) contains the same kind of errors as Eq (5) does. Fortunately, the standard errors are relatively small. Therefore, m1 is the estimated rate of temporal change in pHinsitu and m2βT is the deviation largely due to the various sampling distributions under the strong seasonal pHinsitu variability. In addition, any parameter, having linearity with pHinsitu and having no observable temporal collinearities with t or pHinsitu, can technically be used in the same way as T does to quantify the deviation of βpH due to various sampling distributions under strong seasonal pHinsitu variability.

Linearity between β pH vs. β T and the deviation (m2 β T)

To quantify Eq (6), observations from five time series with varying numbers of observation years (9–23 years) are examined with the average βpH of 31 time series stations along the 137 oE line in the winter and summer (Fig. 4; detailed information are shown in Supplementary Fig. S1 Supplementary Information Table S1). Determined by SLR on the observed pHinsitu and T time series, the results show that the observed βpH values vary from +1.61 × 10−3 pH unit yr−1 to −2.5 × 10−3 pH unit yr−1. The slope (m2) of the regression line of −9.7 ± 0.7 × 10−3 pH unit °C−1 is consistent with the MLR determined m2, with an average of −8.7 ± 1.4 × 10−3 pH unit °C−1 for the five time series (Table 1). Although no reported average values in m1 and m2 are given, the result in Fig. 4 shows that the average βpH values at the 137 oE likely have m1 and m2 values similar to that of the other five time series.

Figure 4
figure 4

The observed surface seawater βpH vs. βT at various time series stations. The solid line shows the linear regression result. The pink solid square shows the expected acidification rate based on the assumption of air-sea CO2 equilibrium. The vertical dashed line represents zero βT. See the text for the names of the time series studies.

Full size image

Assuming an air-sea CO2 equilibrium, the surface ocean has slight regionally distributed differences in acidification rates owing to differences in seawater carbonate chemistry18. Based on thermodynamic calculations using the average carbonate parameters at each station, the acidification rate among the five time series stations and the 137 oE line are similar between −1.6 × 10−3 and −1.8 × 10−3 pH unit yr−1. As βpH is negatively correlated with βT, the deviation (m2βT) is equal to zero when βT is zero. Indeed, the regression line in Fig. 4 yields a βpH of −1.5 ± 0.1 × 10−3 pH unit yr−1 at βT = 0. This value agrees with the MLR-determined m1 of −1.6 ± 0.1 × 10−3 pH unit yr−1 and with the expected acidification rate of −1.6 × 10−3 ~ −1.8 × 10−3 pH unit yr−1, when considering only the observed global atmospheric pCO2 increments during the studied periods. Notably, the assumed, extreme case of uneven sampling distribution shown in Fig. 2 agrees well with the regression line (Fig. 4, red circle). As shown for the first time and having the shortest observation years and the largest deviation from the expected acidification rate (positive βpH), the SEATS site is used as an example to illustrate the influence of sampling distribution on the βpH in the following discussion. At the SEATS site, the observed βpH, m2 and βT of the nine year study are 0.21 × 10−3 pH unit yr−1, −7.7 × 10−3 pH unit °C−1 and −0.23 °C yr−1, respectively (Fig. 4, Table 1). It is important to note that the estimated rate of temporal change in pHinsitu (m1) is βpH – m2βT = (0.21 × 10−3pH unit yr−1)−(−7.7 × 10−3 pH unit oC−1) × (−0.23 °C yr−1) = −1.6 × 10−3 pH unit yr−1, which is statistically indistinguishable from the expected acidification rate, assuming the air-sea CO2 equilibrium.

Discussion

The results above show that the βpH among the studied time-series, at βT = 0, reflects the increasing atmospheric CO2, while deviations from that value are mainly due to different deviations of m2βT under various distributions of sampling points in the studied areas. In short, βpH at ESTOC with only selected data (an extreme case of uneven sampling distribution) and at SEATS and CARIACO, deviates so much from the expected acidification rate due to increasing atmospheric CO2 because of the deviation indicated by the apparent decrease in the observed T at ESTOC and SEATS (positive deviation) and the apparent increase in T at CARIACO (negative deviation). The long-term time-series, such as BATS and HOT, have the βpH statistically indistinguishable from the expected rate of −1.6 × 10−3 ~ −1.8 × 10−3 pH unit yr−1 assuming air-sea CO2 equilibrium. That is, although long-term time series may have evenly distributed datasets and negligible deviations in βpH, shorter time series may suffer large errors which are correctable by this method.

It has been suggested that the anomaly (observed pHinsitu or T minus climatology mean) be used to reduce seasonal variations, as it has been reported that βpH is sensitive to these2. Comparing deseasoned βpH from Bates, et al.2 and deseasoned βT from Astor, et al.9, based on pHinsitu and T anomalies at the CARIACO site, the deseasoned βpH is similar to the observed βpH that contains a deviation of m2βT (Fig. 4, green circle). The reason why the βpH and the deseasoned βpH contain similar deviations needs further investigation.

The observed slope of −9.7 ± 0.7 × 10−3 pH unit °C−1 shown in Fig. 4 can only be due to changing seawater chemistry and not to deviations, if the following two conditions are satisfied. First, βT is a fair estimation of changing T. Second, the change in pHinsitu under long-term change in T is equal to the m2 of −9.7 ± 0.7 × 10−3 pH unit °C−1. However, the second assumption contradicts with our thermodynamic calculation that the change in pHinsitu as T changes (m2) is in fact just approximately −1.1 × 10−3 pH unit oC−1, which is based on the assumptions of air-sea CO2 equilibrium and constant TA. This is consistent with recent studies showing that the change in surface ocean pHinsitu under long-term change in T (m2) is close to zero15,16. Seasonal changes in T are much larger than long-term changes in T, such that the observed m2 of −8.7 ± 1.4 × 10−3 and −9.7 ± 0.7 × 10−3 pH unit °C−1 shown in Table 1 and Fig. 4, respectively, largely represents the values of seasonal change in pHinsitu as T changes. Notably, uneven sampling generates a deviation in βT as well. Although deducing the deviation in βT under uneven sampling distribution is beyond the scope of this study, this does not change a fact that based on Eq (6), in the case of uneven distribution of sampling the slope of βpH vs. βT plot should be approximately equal to m2 = −8.7 ± 1.4 × 10−3pH °C−1. In contrast, the slope of fair estimation is expected to be just about −1.1 × 10−3 pH unit °C−1. The factor of −9.7 ± 0.7 × 10−3 βT shown in Fig. 4 is actually a deviation of the acidification rate largely owing to various distributions of sampling points, rather than an effect on the acidification rate due to warming or cooling.

In the marine carbonate system, the TA concentration should remain constant if the change in pHinsitu is solely caused by the air-sea CO2 exchange. To examine the temporal change in the TA concentration, it is normalized to a particular salinity (S), such as at S = 35 (nTA = TA/S×35), to reduce the influences of precipitation or evaporation on the TA concentration. At the ESTOC and the 137 °E line, the nTA concentrations have reportedly remained unchanged over the study regions and periods8,12. The rates of temporal change in nTA (βnTA) among the rest of the four sites are also found to be insignificant within the uncertainties of the data (see Supplementary Table S1 for detail). Further, changes in TA through precipitation or evaporation have insignificant influences on βpH among those time series stations. For instance, at the HOT site, a change in TA as S changes (11 × 10−3 S unit yr−1) yields a βpH of only 0.091 × 10−3 pH unit yr−1. Such a magnitude of change is less than 5% of the reported βpH at the HOT site. In fact, the constant nTA among the studied time series leads to the conclusion that m1 of −1.6 × 10−3 pH unit yr−1 in Table 1 and −1.5 × 10−3 pH unit yr−1 in Fig. 4 indeed reflect the increasing atmospheric CO2 concentration.

To conclude, long-term time series with even sampling distribution contains negligible deviation. However, a short duration of observation years, low sampling frequency and high seasonal variations in the pHinsitu time series tend to result in a temporal change with a large deviation. For instance, due to its high seasonal variation in pHinsitu and short duration, in the SEATS study, the temporal increase in pHinsitu due to uneven distributions of data is large enough to compensate for decreasing pHinsitu caused by increased atmospheric CO2. T is used as an example to quantify the deviation. Once the deviation is eliminated, seawaters at the studied stations are acidified at a rate statistically indistinguishable from the expected rate under the air-sea CO2 equilibrium. Our results indicate that the differences in βpH among the studied time series are due to deviations using SLR. As our approach is not so limited by the number of observation years or the sampling distribution, time series studies with high seasonal variations and short duration times can benefit from our approach, achieving a fair estimation of the rate of temporal change in pHinsitu.

Methods

The data sets (see Fig. 1 for station locations) used in this study are from five published time series studies (the BATS6, CARIACO2,9, ESTOC8,19, HOT7 and the 137°E12) and an unpublished, open access pH time series study (the SEATS). In this study, the up-to-date CO2 System Calculations Program (version 2.1) developed by Pierrot, et al.20 and the recommended dissociate constants of carbonate chemistry for best practices21,22 were used to calculate the carbonate system. The pH is in the total scale at the in situ temperature. Numbers are expressed as the value ± one standard error.

Additional Information

How to cite this article: Lui, H. K. and Chen, C. T. A. Deducing acidification rates based on short-term time series. Sci. Rep. 5, 11517; doi: 10.1038/srep11517 (2015).