# Applications of Higher-Order Optimal Newton Secant Iterative Methods in Ocean Acidification and Investigation of Long-Run  Implications ofCO2 Emissions on Alkalinity of Seawater

Authors: D. K. R. Babajee; V. C. Jaunky
Journal: ISRN Applied Mathematics
 (2013)

---

## Abstract

The Newton secant method is a third-order iterative nonlinear solver. It requires two function and one first derivative evaluations. However, it is not optimal as it does not satisfy the Kung-Traub conjecture. In this work, we derive an optimal fourth-order Newton secant method with the same number of function evaluations using weight functions and we show that it is a member of the King family of fourth-order methods. We also obtain an eighth-order optimal Newton-secant method. We prove the local convergence of the methods. We apply the methods to solve a fourth-order polynomial arising in ocean acidifications and study their dynamics. We use the data of CO2 available from the National Oceanic and Atmospheric Administration from 1959 to 2012 and calculate the pH of the oceans for these years. Finally we further investigate the long-run implications of CO2 emissions on alkalinity of seawater using fully modified ordinary least squares (FMOLS) and dynamic OLS (DOLS). Our findings reveal that a one-percent increase in CO2 emissions will lead to a reduction in seawater alkalinity of 0.85 percent in the long run.


---

## Body

## 1. Introduction

Recent advancements in the study of higher-order multipoint methods have made this field of research very active. Much literature on the multipoint Newton-like methods for function of one variable and their convergence analysis can be found in [1] and the historical developments of the methods in [2]. Newton secant method is a third-order two-point method and it was rediscovered in [3] as a leapfrog Newton method. However, it is not optimal because the order of an optimal method with 3 function evaluations should be 4 according to the Kung-Traub conjecture. In this work, we derive an optimal fourth-order Newton secant method with same number of function evaluations using weight functions and we show that it is a member of the King family of fourth-order methods. We also obtain an eighth-order optimal Newton secant method. We prove the local convergence of the methods. We apply the methods to solve a fourth-order polynomial arising in ocean acidifications and study their dynamics. We use the data of CO2 available from the National Oceanic and Atmospheric Administration from 1959 to 2012 and calculate the pH of the oceans for these years. Finally, we further investigate the long-run implications of CO2 emissions on alkalinity of seawater using fully modified ordinary least squares (FMOLS) and dynamic OLS (DOLS).


## 2. Developments of the Methods

Letxn+1=ψ(xn) define an iterative function (IF).Definition 1 (see [4]).
If the sequence{xn} tends to a limit x* in such a way that
(1)limn→∞xn+1-x*(xn-x*)p=C
for p≥1, then the order of convergence of the sequence is said to be p, and C is known as the asymptotic error constant. If p=1, p=2, or p=3, then the convergence is said to be linear, quadratic, or cubic, respectively.
Leten=xn-x*, and then the relation
(2)en+1=Cenp+O(enp+1)=O(enp).
is called the error equation. The value of p is called the order of convergence of the method.Definition 2 (see [5]).
The efficiency index is given by(3)EI=p1/d,
where d is the total number of new function evaluations (the values of f and its derivatives) per iteration.Letxn+1 be determined by new information at xn,ϕ1(xn),…,ϕi(xn), i≥1.No old information is reused. Thus,(4)xn+1=ψ(xn,ϕ1(xn),…,ϕi(xn)).
Then ψ is called a multipoint IF without memory.Kung-Traub Conjecture (see [6]). Let ψ be an IF without memory with d evaluations. Then
(5)p(ψ)≤popt=2d-1,
where popt is the maximum order.The Newton (also called Newton-Raphson) IF (2nd NR) is given by(6)ψ2ndNR(x)=x-u(x),u(x)=f(x)f′(x).
The 2nd NR IF is one-point IF with 2 function evaluations and it satisfies the Kung-Traub conjecture d=2.The Halley IF (3rd Hal) is given by(7)ψ3rdHal(x)=x-u(x)1-c2(x)u(x),c2(x)=f′′(x)2f′(x).
It is one-point IF with 3 function evaluations and EI3rdHal=1.44>EI2ndNR=1.41. However, we need to calculate the second derivatives which can be computationally expensive for complex functions. A remedy to this is the Newton-secant IF (3rd NS) which can be written as
(8)ψ3rdNS(x)=x-f(x)f(x)-f[ψ2ndNR(x)]u(x).
It has the same efficiency as the 3rd Hal IF but it requires 2f and 1f′ evaluations, hence no second derivatives, and is a also variant of 3rd Hal IF [7]. However, it does not satisfy the Kung-Traub conjecture. In this work, we develop a fourth-order Newton-secant IF with 3 functions evaluations using weight functions. The 3rd NS's IF can be written as
(9)ψ3rdNS(x)=x-u(x)11-t1,t1=f[ψ2ndNR(x)]f(x).
The 4th NS IF can be given by
(10)ψ4thNS(x)=x-u(x)1+t121-t1.
Since EI4thNS=41/3=1.59, the 4th NS has a higher efficiency index. It is remarkable by just multiplying the last term of (9) by the weight function, 1+t12,  and we could increase the order from 3 to 4. This IF is similar to the member β=1 of the one-parameter King family of fourth-order IF [8] given by
(11)ψ4thFK(x)=ψ2ndNR(x)-f(ψ2ndNR(x))f′(x)1+βt11+(β-2)t1.
Based on King-type family higher-order IFs [9], an eighth-order 8th NS IF can be given by
(12)ψ8thNS(x)=x-1f′(x)(1+t1+1.5t121-t1+0.5t12+4t2+t3),
where
(13)t2=f(ψ4thNS(x))f(x),t3=f(ψ4thNS(x))f(ψ2ndNR(x)).
The 8th NS IF is an optimal IF which satisfies the Kung-Traub conjecture with EI8thNS=51/4=1.68 which is the highest efficiency index among all IFs considered in this work.


## 3. Convergence Analysis

Theorem 3.
Let a sufficiently smooth functionf:D⊂ℝ→ℝ have a simple root x* in the open interval D. Then the 4th NS IF (10) is of local fourth-order convergence and the 8th NS IF (12) is of local eighth-order convergence.Proof.
Letcj=f(j)(x*)/j!f′(x*),  j=2,3,4.
Using the Taylor series and the symbolic software such as, Maple we have(14)f(x)=f′(x*)[en+c2en2+c3en3+c4en4+⋯],(15)f′(x)=f′(x*)[1+2c2en+3c3en2+4c4en3+⋯],
so that
(16)u(x)=en-c2en2+2(c22-c3)en3+(7c2c3-4c23-3c4)en4+⋯,(17)ψ2ndNR(x)-x*=c2en2-2(c22-c3)en3-(7c2c3-4c23-3c4)en4+⋯.
Now, the Taylor expansion of f(y) about x* gives
(18)f(y)=f′(x*)×[(y-x*)+c2(y-x*)2+c3(y-x*)3+c4(y-x*)4+⋯].
Using (14), (18), and (17), we have
(19)t1=f[ψ2ndNR(x)]f(x)=c2en+(2c3-3c22)en2+(3c4-10c2c3+8c23)en3+(-14c2c4+37c3c22-20c24-8c32+4c5)en4+⋯,
so that by using computer algebra software such as Maple we get
(20)ψ3thNS(x)-x*=c22en3+⋯,ψ4thNS(x)-x*=c2(-c3+3c22)en4+⋯.
Similarly, we have
(21)t2=c2(-c3+3c22)en3+(-2c2c4-2c32+21c3c22-21c24)en4+⋯,t3=(-c3+3c22)en2+(-2c4+12c2c3-12c23)en3+⋯,
so that finally we get
(22)ψ8thNS(x)-x*=(12c2(-c3+3c22)(65c24-34c3c22+2c2c4+2c32))en8+⋯.


## 4. Ocean Acidification

### 4.1. Introduction [2]

The accumulation of greenhouse gases (GHGs) in the Earth's atmosphere is now a major topic of discussion to anticipate changes in the Earth's climate. The GHGs cause a reduction in the reradiation of energy from the Sun back into the outer space. Since less energy leaves the Earth's atmosphere, heating of the atmosphere results as a manifest in a temperature rise [12]. This temperature rise, the so-called global warming, is in turn a driving force for climate change. CO2 is the major GHG, with increasing levels primarily from the burning of fossil fuels. Thus, changes in the CO2 level or concentration in the Earth's atmosphere are of paramount importance in understanding anticipated warming and climate change. A second aspect of CO2 accumulation in the atmosphere that is not as generally recognized and appreciated as temperature rise is the accumulation of carbon (from CO2) in the oceans that leads to ocean acidification. CO2 dissolves in ocean water and undergoes a series of chemical changes that ultimately leads to increased hydrogen ion concentration, denoted subsequently as [H+], and thus acidification (see [12, 13]). This increase in [H+] is manifest as a decrease in pH; note that [H+] and pH move in opposite directions due to the following basic relation:
(23)pH=-log10[H+].
Ocean acidification is the name given to the ongoing decrease in the pH of the Earth's oceans, caused by their uptake of anthropogenic carbon dioxide from the atmosphere [14]. Between 1751 and 1994, surface ocean pH is estimated to have decreased from approximately 8.179 to 8.104 (a change of −0.075) [15]. A decrease in ocean pH of 0.1 units corresponds to a 30%  increase in the concentration of H+ in seawater, assuming that alkalinity and temperature remain constant [16, p. 406]. There is about fifty times as much carbon dissolved in the oceans in the form of CO2 and carbonic acid, bicarbonate, and carbonate ions as that in the atmosphere. The oceans act as an enormous carbon sink and have taken up about a third of CO2 emitted by human activity [17]. Most of the CO2 taken up by the ocean forms carbonic acid in equilibrium with bicarbonate and carbonate ions. Some is consumed in photosynthesis by organisms in the water, and a small proportion of that sinks and leaves the carbon cycle. Increased CO2 in the atmosphere has led to decreasing alkalinity of seawater and there is concern that this may adversely affect organisms living in the water. In particular, with decreasing alkalinity, the availability of carbonates for forming shells decreases [18, p. 125], although there is evidence of increased shell production by certain species under increased CO2 content [19].


### 4.2. Ocean Chemistry [2]

We begin withCO2 dissolving in H2O to form carbonic acid, H2CO3 as
(24)CO2+H2O⇌H2CO3.
The double arrow ⇌ denotes a reversible chemical reaction (a reaction that can proceed either forward to produce H2CO3 or backward to produce CO2 and H2O).A common convention is to take[CO2] as the dissolved CO2, denoted as [CO2]aq, plus the carbonic acid [H2CO3]aq.Carbonic acid is a weak acid which in turn dissociate into bicarbonate ions,HCO3-, as
(25)H2CO3⇌H++HCO3-.
Bicarbonate ions in turn dissociates into carbonate ions, CO32-(26)HCO3-⇌H++CO32-.
The reactions (25) and (26) produce hydrogen ions and therefore contribute to acidification.H
2
O also dissociates to produce hydrogen ions as
(27)H2O⇌H++OH-.
Additionally, boron hydroxide in seawater dissociates to produce hydrogen ions as
(28)B(OH)3+H2O⇌H++B(OH)4-.
In this work, we do not consider other compounds in the oceans that also dissociate to produce hydrogen ions. We explain how to compute [H+] and pH from the reactions (24) to (28). For the following analysis, we use the equilibrium constants of Bacastow and Keeling [20] expressed in the units of moles/litre. The relation between gaseous and liquid CO2 is
(29)K0=[CO2]Pt=3.347(-5),
where [CO2] is the sum of the dissolved CO2 and carbonic acid and Pt is the gas phase CO2 partial pressure in ppm measured by the National Oceanic and Atmospheric Administration (NOAA) at the Mauna Loa Observatory, Hawaii [21], and a(-b) denotes a×10-b.For reaction (25),
(30)K1=[H+][HCO3-][CO2]=9.747(-7).For reaction (26),
(31)K2=[H+][CO32-][HCO3-]=8.501(-10).For reaction (27),
(32)KW=[H+][OH-]=6.46(-15).For reaction (28),
(33)KB=[H+][B(OH)4-][B(OH)3]=1.881(-9).
The alkalinity, A, which expresses the electrical neutrality of ocean water is defined as
(34)A=[HCO3-]+2[CO32-]+[B(OH)4-]+[OH-]-[H+].
We can assume that the values of A do not change with time [20]. From (29), we have
(35)[CO2]=K0Pt.
From (30) and (35), we have
(36)[HCO3-]=K1[CO2][H+]=K0K1Pt[H+].
Similarly, we obtain
(37)[CO32-]=K2[HCO3-][H+]=K0K1K2Pt[H+]2
from (31) and (36).Using(38)B=[B(OH)3]+[B(OH)4-]
in (33), we have
(39)[B(OH)4-]=BKBKB+[H+].
Substituting (32) and (35)–(39) into (34), we have
(40)A=K0K1Pt[H+]+2K0K1K2Pt[H+]2+BKBKB+[H+]+KW[H+]-[H+]
which simplifies to the solution of a fourth-order polynomial given by
(41)p([H+])=∑n=04‍Dn[H+]n,
where
(42)D0=2K0K1K2PtKB,D1=K0K1PtKB+2K0K1K2Pt+KWKB,D2=K0K1Pt+BKB+KW-AKB,D3=-KB-A,D4=-1.
We use A=2.050 [22, p. 334] and B=0.409 [20, p. 131].


### 4.3. Dynamic Behaviour

We next study the dynamic behaviour of the methods in the complex planeℤ to find the best starting points. For a given value of Pt, polynomial p([H+]) in (41) has one positive real root (the one we are seeking), one negative real root, and two complex roots. Since these solutions have very small values except the negative one, it is difficult to study their polynomiography [23]. Instead, we consider the change of variable
(43)z=1[H+],z∈ℤ,
and then the pH=log10z. We require to find the positive real solution of another fourth-order polynomial:
(44)p(z)=∑n=04‍D4-nzn.We draw the polynomiographs ofp([H+]). Let z0=x+iy, and let x,y∈ℝ be the initial point. A square grid of 65536 points, composed of 256 columns and 256 rows corresponding to the pixels of a computer display, would represent a region of the complex plane [24]. We consider the square ℝ×ℝ=[-1(9),1(9)]×[-1(9),1(9)]. Each grid point is used as a starting value z0 of the sequence zk+1=ψIF(zk) and the number of iterations until convergence is counted for each grid point. We assign different colours to each root zj*, j=1,2,3,4 of p([H+]) if |zj*-zk|<1(-4), in at most 25 iterations In this way, the basin of attraction for each root would be assigned a characteristic colour. The common boundaries of these basins of attraction constitute the Julia set of the IF If the iterates do not satisfy the above criterion for convergence, we assign the dark blue colour.Figures1 and 2 show the polynomiographs of the 2nd NR, 3rd NS, 4th NS and 8th NS methods, respectively. In this case, the positive root of the polynomial, z*=1.262801212073384(8) (coloured brownish yellow), corresponds to the solution [H+]*=7.918902757133942(-9). It can be shown that there are diverging points for the higher-order Newton-secant methods and that the 2nd NR method has the largest basins of attraction for the positive root among the 4 methods. But we are using the dynamics of the methods to find a suitable starting point for the higher-order Newton-secant methods so that we can make use of their higher-order convergence. Figures 3 and 4 show the basins of attractions on the real line of the 2nd NR, 3rd NS, 4th NS, and 8th NS methods, respectively, for the positive root of p(z). They reveal that the 2nd NR will converge for a starting point z0>0.7(8). As the order of the method increase, the basins of attraction decrease and higher-order Newton-secant methods have difficulty to converge for some starting points. We also find that all methods will converge for the starting point z0=1.25(8) or [H+]0=8(-9).Polynomiographs of the 2nd NR and 3rd NS methods for the polynomialp(z) with Pt=393.81.
(a)
2nd NR
(b)
3rd NSPolynomiographs of the 4th NS and 8th NS methods for the polynomialp(z) with Pt=393.81.
(a)
4th NS
(b)
8th NSBasins of attractions on the real line of the 2nd NR and 3rd NS methods for the positive root  ofp(z).
(a)
2nd NR
(b)
3rd NSBasins of attractions on the real line of the 4th NS and 8th NS methods for the positive root ofp(z).
(a)
4th NS
(b)
8th NS


### 4.4. Numerical Experiments and Results

We use the data available from NOOA to calculate the pH of the ocean from 1959 to 2012. We use a common starting point[H]0+=8(-9) for each Pt and stop the methods whenever |[H]k+1+-[H]k+|<1(-12) in at most 25 iterations. The approximate solutions are calculated correctly to 16 digits in MATLAB. We denote by Ns the number of successful points and by μ as the mean iteration number for the converging points. Table 2 gives a comparison in which we observe that the 3 methods successfully converge to the required root but the 8th NS method has a few diverging points. The 4th NS method is the most effective with the lowest mean iteration number and all converging points. Table 1 shows the calculated pH from 1959 to 2012. Figure 5 shows the variation of CO2 and pH with time. We observe that as the CO2 increases, the pH decreases.Table 1
pH of oceans using thePt from NOAA from 1959 to 2012.
Time
P
t
pH
Time
P
t
pH
1959
315.98
8.1794
1986
347.19
8.1463
1960
316.91
8.1784
1987
348.98
8.1444
1961
317.64
8.1776
1988
351.45
8.1419
1962
318.45
8.1767
1989
352.90
8.1405
1963
318.99
8.1761
1990
354.16
8.1392
1964
319.62
8.1754
1991
355.48
8.1379
1965
320.04
8.1749
1992
356.27
8.1371
1966
321.38
8.1735
1993
356.95
8.1364
1967
322.16
8.1726
1994
358.64
8.1347
1968
323.04
8.1717
1995
360.62
8.1328
1969
324.62
8.1699
1996
362.36
8.1311
1970
325.68
8.1688
1997
363.47
8.1300
1971
326.32
8.1681
1998
366.50
8.1270
1972
327.45
8.1669
1999
368.14
8.1254
1973
329.68
8.1645
2000
369.40
8.1242
1974
330.17
8.1640
2001
371.07
8.1226
1975
331.08
8.1630
2002
373.17
8.1206
1976
332.05
8.1620
2003
375.78
8.1181
1977
333.78
8.1602
2004
377.52
8.1165
1978
335.41
8.1584
2005
379.76
8.1144
1979
336.78
8.1570
2006
381.85
8.1124
1980
338.68
8.1550
2007
383.71
8.1107
1981
341.11
8.1525
2008
385.57
8.1089
1982
341.22
8.1524
2009
387.35
8.1073
1983
342.84
8.1507
2010
389.85
8.1050
1984
344.41
8.1491
2011
391.62
8.1033
1985
345.87
8.1476
2012
393.81
8.1013Table 2
Comparison of successful starting point and mean iteration number for each method.
Method
N
s
μ
2nd NR
54
3.7037
3th NS
54
2.8704
4th NS
54
2.7778
8th NS
47
2.1852Figure 5
Variation of CO2 and pH with time.


### 4.5. Empirical Analysis of Impact of CO2 on Alkalinity of Seawater

To empirically test the impact ofCO2 in the atmosphere on the alkalinity of seawater, we set up the following generalized equation:
(45)pH=f(CO2,ϵ),
where ϵ is the error term. The concept of cointegration as per Engle and Granger [25] is used to investigate any long-run relationship between nonstationary variables. Time-series data such as pH and CO2 tend to be nonstationary in levels. If a series is stationary, then the probability laws controlling its process are stable over time, that is, in statistical equilibrium [26]. In contrast, series having a unit root are nonstationary. Shocks have a unit root and can, in part, change the long-run level of the time series permanently. Per se, a series is said to be integrated of order υ or I(υ) if it were to be different by υ times to become stationary. A stationary process is a series which follows an I(0) process. To run the model, the logarithm of base 10 of the variables is taken. As a prerequisite of the cointegration test, the unit root properties of the two series are investigated. The augmented Dickey-Fuller (ADF) test as proposed by Dickey and Fuller [27] and the DF-GLS test as per Elliott et al. [28] for the null of a unit root are considered. The DF-GLS test is a modified ADF test and tends to be a more asymptotically powerful test. These tests apply regressions which include a constant term only, while the other contain both a constant term and a time trend. Time series data tend to exhibit a trend over time and hence it is more appropriate to consider a regression with both a constant term and a trend. In contrast, first differencing is likely to remove any deterministic trends. Hence, the regression should include a constant only. In general, time-series data tends to be nonstationary and I(1). Both series must be integrated of the same order to validate a cointegrating relationship. The Johansen cointegration test [29] is conducted within a vector autoregression (VAR) structure and it involves two log-likelihood ratio (LR) test statistics, namely, the maximum eigenvalue (λ-max) and trace (Tr) statistics. Once a cointegrating relationship is established, long-run estimates can be computed via the fully modified ordinary least squares (FMOLS) and dynamic OLS (DOLS) of Phillips and Hansen [30] and Stock and Watson [31], respectively. Table 3 shows the results of the unit root tests. Both series are found to be nonstationary. The ADF test statistics illustrate an I(1) process for both series only when a trend is considered in the testing framework. However, when testing for a unit root using first-differenced data, the trend should be excluded. The DF-GLS confirms our a priori expectation. Both series are found to be I(1) for both deterministics. Table 4 reports the cointegration test statistics. According to the null hypothesis for the λ-max and Tr tests, there are at most r cointegrating vectors, whereas the alternative hypotheses are r+1 and at least r+1 for the λ-max and Tr statistics, respectively. As per the λ-max statistics, the null hypothesis of r=0 is rejected in favour of r=1. A similar result is found when referring to the Tr statistics as the null hypothesis of r=0 is rejected in favour of r≥1. The computed test statistics are 23.75 and 28.53 for the λ-max and Tr tests, respectively. The null hypothesis of no cointegration is rejected at 5% level. Furthermore, the null hypothesis of at most one cointegrating vector (r≤1) is in no case rejected in both cases. In sum, these findings provide evidence of a long-run equilibrium relationship between pH and CO2. Given the presence of a cointegrating vector, the long-run elasticity can now be computed and is reported in Table 5. The FMOLS and DOLS methods are robust single equation approaches which can correct for endogeneity bias and serial correlation (The computed test statistic for serial correlation according to Durbin and Watson [32] is d-statistic (2,54)=0.021. This reveals positive serial correlation) in a semiparametric and parametric way, respectively. CO2 in the atmosphere has a statistically significant negative impact on the alkalinity of seawater and the long-run elasticities from both methods tend to coincide. For instance, a one-percent increase in CO2 emissions will generate to a reduction in seawater alkalinity of 0.85 percent in the long run.Table 3
Unit root tests.
Series
ADF
DF-GLS
With constant and without trend
With constant and with trend
With constant and without trend
With constant and with trend
log
10
pH
3.412
[
0
]
-
1.800
[
0
]
-
0.121
[
4
]
-
0.147
[
1
]
log
10
C
O
2
3.266
[
0
]
-
1.894
[
0
]
-
0.130
[
4
]
-
0.188
[
1
]
Δ
log
10
pH
-
2.291
[
3
]
-
6.104
[
1
]
*
-
2.352
[
2
]
*
*
-
4.486
[
1
]
*
Δ
log
10
C
O
2
-
2.347
[
3
]
-
6.101
[
1
]
*
-
2.423
[
2
]
*
*
-
6.241
[
1
]
*
Note: to select the order of lag, we start with a maximum lag length of 4 and pare it down as per the Akaike information criterion (AIC). There is no general rule on how to choose the maximum lag to start with. The bandwidth and maximum lag length are chosen according to the Bartlett kernel which is equal to4(T/100)2/9≈4, where T=54. The optimal lag length is given in square brackets. The MacKinnon critical values [10] for the ADF unit root tests with a constant and without a time are -3.59, -2.94, and -2.60 at 1%, 5%, and 10% significance level,respectively, while those with a constant and a time trend are -4.17, -3.51, and -3.19, respectively. DF-GLS critical values without trend at 1%, 5%, and 10% levels are -2.62, -2.26, and -1.95 and with a trend are -3.76, -3.17, and -2.87, respectively. The optimal lag is chosen according to the Akaike information criterion (AIC) and Schwarz Bayesian criterion for the ADF and DF-GLS tests, respectively.                   * and ** denote 1% and 5% significance level correspondingly.Table 4
Johansen cointegration test.
LR test
Hypothesis
Statistics
95
% critical values
90
% critical values
Null
Alternative
λ-max
r
=
0
r
=
1
23.75
2
*
*
18.330
16.280
r
≤
1
r
=
2
4.773
11.540
9.750
Tr
r
=
0
r
≥
1
28.52
5
*
*
23.830
21.230
r
≤
1
r
=
2
4.773
11.540
9.750
Note: the test is conducted with unrestricted constants and trends in the VAR model.r is the number of cointegrating vectors. The optimal lag length is set to  4 according to the AIC.Table 5
Long-run estimators.
Series
Dependentlog10pH
FMOLS
DOLS
Coefficient
Standard deviation
Coefficient
Standard deviation
log
10
C
O
2
-
0.84
5
*
0.003
-
0.84
9
*
0.009
Note: a constant and time trend are  included in each model. The critical values of the two-tailedt-statistics test at 1%, 5%, and 10% significance levels are 2.326, 1.645, and 1.282, respectively. The maximum lag/lead is set to 2 [11].


## 4.1. Introduction [2]

The accumulation of greenhouse gases (GHGs) in the Earth's atmosphere is now a major topic of discussion to anticipate changes in the Earth's climate. The GHGs cause a reduction in the reradiation of energy from the Sun back into the outer space. Since less energy leaves the Earth's atmosphere, heating of the atmosphere results as a manifest in a temperature rise [12]. This temperature rise, the so-called global warming, is in turn a driving force for climate change. CO2 is the major GHG, with increasing levels primarily from the burning of fossil fuels. Thus, changes in the CO2 level or concentration in the Earth's atmosphere are of paramount importance in understanding anticipated warming and climate change. A second aspect of CO2 accumulation in the atmosphere that is not as generally recognized and appreciated as temperature rise is the accumulation of carbon (from CO2) in the oceans that leads to ocean acidification. CO2 dissolves in ocean water and undergoes a series of chemical changes that ultimately leads to increased hydrogen ion concentration, denoted subsequently as [H+], and thus acidification (see [12, 13]). This increase in [H+] is manifest as a decrease in pH; note that [H+] and pH move in opposite directions due to the following basic relation:
(23)pH=-log10[H+].
Ocean acidification is the name given to the ongoing decrease in the pH of the Earth's oceans, caused by their uptake of anthropogenic carbon dioxide from the atmosphere [14]. Between 1751 and 1994, surface ocean pH is estimated to have decreased from approximately 8.179 to 8.104 (a change of −0.075) [15]. A decrease in ocean pH of 0.1 units corresponds to a 30%  increase in the concentration of H+ in seawater, assuming that alkalinity and temperature remain constant [16, p. 406]. There is about fifty times as much carbon dissolved in the oceans in the form of CO2 and carbonic acid, bicarbonate, and carbonate ions as that in the atmosphere. The oceans act as an enormous carbon sink and have taken up about a third of CO2 emitted by human activity [17]. Most of the CO2 taken up by the ocean forms carbonic acid in equilibrium with bicarbonate and carbonate ions. Some is consumed in photosynthesis by organisms in the water, and a small proportion of that sinks and leaves the carbon cycle. Increased CO2 in the atmosphere has led to decreasing alkalinity of seawater and there is concern that this may adversely affect organisms living in the water. In particular, with decreasing alkalinity, the availability of carbonates for forming shells decreases [18, p. 125], although there is evidence of increased shell production by certain species under increased CO2 content [19].


## 4.2. Ocean Chemistry [2]

We begin withCO2 dissolving in H2O to form carbonic acid, H2CO3 as
(24)CO2+H2O⇌H2CO3.
The double arrow ⇌ denotes a reversible chemical reaction (a reaction that can proceed either forward to produce H2CO3 or backward to produce CO2 and H2O).A common convention is to take[CO2] as the dissolved CO2, denoted as [CO2]aq, plus the carbonic acid [H2CO3]aq.Carbonic acid is a weak acid which in turn dissociate into bicarbonate ions,HCO3-, as
(25)H2CO3⇌H++HCO3-.
Bicarbonate ions in turn dissociates into carbonate ions, CO32-(26)HCO3-⇌H++CO32-.
The reactions (25) and (26) produce hydrogen ions and therefore contribute to acidification.H
2
O also dissociates to produce hydrogen ions as
(27)H2O⇌H++OH-.
Additionally, boron hydroxide in seawater dissociates to produce hydrogen ions as
(28)B(OH)3+H2O⇌H++B(OH)4-.
In this work, we do not consider other compounds in the oceans that also dissociate to produce hydrogen ions. We explain how to compute [H+] and pH from the reactions (24) to (28). For the following analysis, we use the equilibrium constants of Bacastow and Keeling [20] expressed in the units of moles/litre. The relation between gaseous and liquid CO2 is
(29)K0=[CO2]Pt=3.347(-5),
where [CO2] is the sum of the dissolved CO2 and carbonic acid and Pt is the gas phase CO2 partial pressure in ppm measured by the National Oceanic and Atmospheric Administration (NOAA) at the Mauna Loa Observatory, Hawaii [21], and a(-b) denotes a×10-b.For reaction (25),
(30)K1=[H+][HCO3-][CO2]=9.747(-7).For reaction (26),
(31)K2=[H+][CO32-][HCO3-]=8.501(-10).For reaction (27),
(32)KW=[H+][OH-]=6.46(-15).For reaction (28),
(33)KB=[H+][B(OH)4-][B(OH)3]=1.881(-9).
The alkalinity, A, which expresses the electrical neutrality of ocean water is defined as
(34)A=[HCO3-]+2[CO32-]+[B(OH)4-]+[OH-]-[H+].
We can assume that the values of A do not change with time [20]. From (29), we have
(35)[CO2]=K0Pt.
From (30) and (35), we have
(36)[HCO3-]=K1[CO2][H+]=K0K1Pt[H+].
Similarly, we obtain
(37)[CO32-]=K2[HCO3-][H+]=K0K1K2Pt[H+]2
from (31) and (36).Using(38)B=[B(OH)3]+[B(OH)4-]
in (33), we have
(39)[B(OH)4-]=BKBKB+[H+].
Substituting (32) and (35)–(39) into (34), we have
(40)A=K0K1Pt[H+]+2K0K1K2Pt[H+]2+BKBKB+[H+]+KW[H+]-[H+]
which simplifies to the solution of a fourth-order polynomial given by
(41)p([H+])=∑n=04‍Dn[H+]n,
where
(42)D0=2K0K1K2PtKB,D1=K0K1PtKB+2K0K1K2Pt+KWKB,D2=K0K1Pt+BKB+KW-AKB,D3=-KB-A,D4=-1.
We use A=2.050 [22, p. 334] and B=0.409 [20, p. 131].


## 4.3. Dynamic Behaviour

We next study the dynamic behaviour of the methods in the complex planeℤ to find the best starting points. For a given value of Pt, polynomial p([H+]) in (41) has one positive real root (the one we are seeking), one negative real root, and two complex roots. Since these solutions have very small values except the negative one, it is difficult to study their polynomiography [23]. Instead, we consider the change of variable
(43)z=1[H+],z∈ℤ,
and then the pH=log10z. We require to find the positive real solution of another fourth-order polynomial:
(44)p(z)=∑n=04‍D4-nzn.We draw the polynomiographs ofp([H+]). Let z0=x+iy, and let x,y∈ℝ be the initial point. A square grid of 65536 points, composed of 256 columns and 256 rows corresponding to the pixels of a computer display, would represent a region of the complex plane [24]. We consider the square ℝ×ℝ=[-1(9),1(9)]×[-1(9),1(9)]. Each grid point is used as a starting value z0 of the sequence zk+1=ψIF(zk) and the number of iterations until convergence is counted for each grid point. We assign different colours to each root zj*, j=1,2,3,4 of p([H+]) if |zj*-zk|<1(-4), in at most 25 iterations In this way, the basin of attraction for each root would be assigned a characteristic colour. The common boundaries of these basins of attraction constitute the Julia set of the IF If the iterates do not satisfy the above criterion for convergence, we assign the dark blue colour.Figures1 and 2 show the polynomiographs of the 2nd NR, 3rd NS, 4th NS and 8th NS methods, respectively. In this case, the positive root of the polynomial, z*=1.262801212073384(8) (coloured brownish yellow), corresponds to the solution [H+]*=7.918902757133942(-9). It can be shown that there are diverging points for the higher-order Newton-secant methods and that the 2nd NR method has the largest basins of attraction for the positive root among the 4 methods. But we are using the dynamics of the methods to find a suitable starting point for the higher-order Newton-secant methods so that we can make use of their higher-order convergence. Figures 3 and 4 show the basins of attractions on the real line of the 2nd NR, 3rd NS, 4th NS, and 8th NS methods, respectively, for the positive root of p(z). They reveal that the 2nd NR will converge for a starting point z0>0.7(8). As the order of the method increase, the basins of attraction decrease and higher-order Newton-secant methods have difficulty to converge for some starting points. We also find that all methods will converge for the starting point z0=1.25(8) or [H+]0=8(-9).Polynomiographs of the 2nd NR and 3rd NS methods for the polynomialp(z) with Pt=393.81.
(a)
2nd NR
(b)
3rd NSPolynomiographs of the 4th NS and 8th NS methods for the polynomialp(z) with Pt=393.81.
(a)
4th NS
(b)
8th NSBasins of attractions on the real line of the 2nd NR and 3rd NS methods for the positive root  ofp(z).
(a)
2nd NR
(b)
3rd NSBasins of attractions on the real line of the 4th NS and 8th NS methods for the positive root ofp(z).
(a)
4th NS
(b)
8th NS


## 4.4. Numerical Experiments and Results

We use the data available from NOOA to calculate the pH of the ocean from 1959 to 2012. We use a common starting point[H]0+=8(-9) for each Pt and stop the methods whenever |[H]k+1+-[H]k+|<1(-12) in at most 25 iterations. The approximate solutions are calculated correctly to 16 digits in MATLAB. We denote by Ns the number of successful points and by μ as the mean iteration number for the converging points. Table 2 gives a comparison in which we observe that the 3 methods successfully converge to the required root but the 8th NS method has a few diverging points. The 4th NS method is the most effective with the lowest mean iteration number and all converging points. Table 1 shows the calculated pH from 1959 to 2012. Figure 5 shows the variation of CO2 and pH with time. We observe that as the CO2 increases, the pH decreases.Table 1
pH of oceans using thePt from NOAA from 1959 to 2012.
Time
P
t
pH
Time
P
t
pH
1959
315.98
8.1794
1986
347.19
8.1463
1960
316.91
8.1784
1987
348.98
8.1444
1961
317.64
8.1776
1988
351.45
8.1419
1962
318.45
8.1767
1989
352.90
8.1405
1963
318.99
8.1761
1990
354.16
8.1392
1964
319.62
8.1754
1991
355.48
8.1379
1965
320.04
8.1749
1992
356.27
8.1371
1966
321.38
8.1735
1993
356.95
8.1364
1967
322.16
8.1726
1994
358.64
8.1347
1968
323.04
8.1717
1995
360.62
8.1328
1969
324.62
8.1699
1996
362.36
8.1311
1970
325.68
8.1688
1997
363.47
8.1300
1971
326.32
8.1681
1998
366.50
8.1270
1972
327.45
8.1669
1999
368.14
8.1254
1973
329.68
8.1645
2000
369.40
8.1242
1974
330.17
8.1640
2001
371.07
8.1226
1975
331.08
8.1630
2002
373.17
8.1206
1976
332.05
8.1620
2003
375.78
8.1181
1977
333.78
8.1602
2004
377.52
8.1165
1978
335.41
8.1584
2005
379.76
8.1144
1979
336.78
8.1570
2006
381.85
8.1124
1980
338.68
8.1550
2007
383.71
8.1107
1981
341.11
8.1525
2008
385.57
8.1089
1982
341.22
8.1524
2009
387.35
8.1073
1983
342.84
8.1507
2010
389.85
8.1050
1984
344.41
8.1491
2011
391.62
8.1033
1985
345.87
8.1476
2012
393.81
8.1013Table 2
Comparison of successful starting point and mean iteration number for each method.
Method
N
s
μ
2nd NR
54
3.7037
3th NS
54
2.8704
4th NS
54
2.7778
8th NS
47
2.1852Figure 5
Variation of CO2 and pH with time.


## 4.5. Empirical Analysis of Impact of CO2 on Alkalinity of Seawater

To empirically test the impact ofCO2 in the atmosphere on the alkalinity of seawater, we set up the following generalized equation:
(45)pH=f(CO2,ϵ),
where ϵ is the error term. The concept of cointegration as per Engle and Granger [25] is used to investigate any long-run relationship between nonstationary variables. Time-series data such as pH and CO2 tend to be nonstationary in levels. If a series is stationary, then the probability laws controlling its process are stable over time, that is, in statistical equilibrium [26]. In contrast, series having a unit root are nonstationary. Shocks have a unit root and can, in part, change the long-run level of the time series permanently. Per se, a series is said to be integrated of order υ or I(υ) if it were to be different by υ times to become stationary. A stationary process is a series which follows an I(0) process. To run the model, the logarithm of base 10 of the variables is taken. As a prerequisite of the cointegration test, the unit root properties of the two series are investigated. The augmented Dickey-Fuller (ADF) test as proposed by Dickey and Fuller [27] and the DF-GLS test as per Elliott et al. [28] for the null of a unit root are considered. The DF-GLS test is a modified ADF test and tends to be a more asymptotically powerful test. These tests apply regressions which include a constant term only, while the other contain both a constant term and a time trend. Time series data tend to exhibit a trend over time and hence it is more appropriate to consider a regression with both a constant term and a trend. In contrast, first differencing is likely to remove any deterministic trends. Hence, the regression should include a constant only. In general, time-series data tends to be nonstationary and I(1). Both series must be integrated of the same order to validate a cointegrating relationship. The Johansen cointegration test [29] is conducted within a vector autoregression (VAR) structure and it involves two log-likelihood ratio (LR) test statistics, namely, the maximum eigenvalue (λ-max) and trace (Tr) statistics. Once a cointegrating relationship is established, long-run estimates can be computed via the fully modified ordinary least squares (FMOLS) and dynamic OLS (DOLS) of Phillips and Hansen [30] and Stock and Watson [31], respectively. Table 3 shows the results of the unit root tests. Both series are found to be nonstationary. The ADF test statistics illustrate an I(1) process for both series only when a trend is considered in the testing framework. However, when testing for a unit root using first-differenced data, the trend should be excluded. The DF-GLS confirms our a priori expectation. Both series are found to be I(1) for both deterministics. Table 4 reports the cointegration test statistics. According to the null hypothesis for the λ-max and Tr tests, there are at most r cointegrating vectors, whereas the alternative hypotheses are r+1 and at least r+1 for the λ-max and Tr statistics, respectively. As per the λ-max statistics, the null hypothesis of r=0 is rejected in favour of r=1. A similar result is found when referring to the Tr statistics as the null hypothesis of r=0 is rejected in favour of r≥1. The computed test statistics are 23.75 and 28.53 for the λ-max and Tr tests, respectively. The null hypothesis of no cointegration is rejected at 5% level. Furthermore, the null hypothesis of at most one cointegrating vector (r≤1) is in no case rejected in both cases. In sum, these findings provide evidence of a long-run equilibrium relationship between pH and CO2. Given the presence of a cointegrating vector, the long-run elasticity can now be computed and is reported in Table 5. The FMOLS and DOLS methods are robust single equation approaches which can correct for endogeneity bias and serial correlation (The computed test statistic for serial correlation according to Durbin and Watson [32] is d-statistic (2,54)=0.021. This reveals positive serial correlation) in a semiparametric and parametric way, respectively. CO2 in the atmosphere has a statistically significant negative impact on the alkalinity of seawater and the long-run elasticities from both methods tend to coincide. For instance, a one-percent increase in CO2 emissions will generate to a reduction in seawater alkalinity of 0.85 percent in the long run.Table 3
Unit root tests.
Series
ADF
DF-GLS
With constant and without trend
With constant and with trend
With constant and without trend
With constant and with trend
log
10
pH
3.412
[
0
]
-
1.800
[
0
]
-
0.121
[
4
]
-
0.147
[
1
]
log
10
C
O
2
3.266
[
0
]
-
1.894
[
0
]
-
0.130
[
4
]
-
0.188
[
1
]
Δ
log
10
pH
-
2.291
[
3
]
-
6.104
[
1
]
*
-
2.352
[
2
]
*
*
-
4.486
[
1
]
*
Δ
log
10
C
O
2
-
2.347
[
3
]
-
6.101
[
1
]
*
-
2.423
[
2
]
*
*
-
6.241
[
1
]
*
Note: to select the order of lag, we start with a maximum lag length of 4 and pare it down as per the Akaike information criterion (AIC). There is no general rule on how to choose the maximum lag to start with. The bandwidth and maximum lag length are chosen according to the Bartlett kernel which is equal to4(T/100)2/9≈4, where T=54. The optimal lag length is given in square brackets. The MacKinnon critical values [10] for the ADF unit root tests with a constant and without a time are -3.59, -2.94, and -2.60 at 1%, 5%, and 10% significance level,respectively, while those with a constant and a time trend are -4.17, -3.51, and -3.19, respectively. DF-GLS critical values without trend at 1%, 5%, and 10% levels are -2.62, -2.26, and -1.95 and with a trend are -3.76, -3.17, and -2.87, respectively. The optimal lag is chosen according to the Akaike information criterion (AIC) and Schwarz Bayesian criterion for the ADF and DF-GLS tests, respectively.                   * and ** denote 1% and 5% significance level correspondingly.Table 4
Johansen cointegration test.
LR test
Hypothesis
Statistics
95
% critical values
90
% critical values
Null
Alternative
λ-max
r
=
0
r
=
1
23.75
2
*
*
18.330
16.280
r
≤
1
r
=
2
4.773
11.540
9.750
Tr
r
=
0
r
≥
1
28.52
5
*
*
23.830
21.230
r
≤
1
r
=
2
4.773
11.540
9.750
Note: the test is conducted with unrestricted constants and trends in the VAR model.r is the number of cointegrating vectors. The optimal lag length is set to  4 according to the AIC.Table 5
Long-run estimators.
Series
Dependentlog10pH
FMOLS
DOLS
Coefficient
Standard deviation
Coefficient
Standard deviation
log
10
C
O
2
-
0.84
5
*
0.003
-
0.84
9
*
0.009
Note: a constant and time trend are  included in each model. The critical values of the two-tailedt-statistics test at 1%, 5%, and 10% significance levels are 2.326, 1.645, and 1.282, respectively. The maximum lag/lead is set to 2 [11].


## 5. Conclusion

We develop an optimal fourth- and eighth-order Newton-secant methods. We study their dynamics in a fourth-order polynomial arising in ocean acidification. We also perform an investigation on the long-run implications ofCO2 emissions on alkalinity of seawater using fully modified ordinary least squares (FMOLS) and dynamic OLS (DOLS). We find that a one-percent increase in CO2 emissions will lead to a reduction in seawater alkalinity of 0.85 percent in the long run. Put differently, a fall in CO2 emissions will lead to an improvement of the quality of seawater and therefore to the sustainability of the marine ecosystem.



---

*Source: 785287-2013-05-29.xml*