Natural variability and anthropogenic effects in a Central Mediterranean core

. We evaluate the contribution of natural variability to the modern decrease in foraminiferal δ 18 O by relying on a 2200-yr-long, high-resolution record of oxygen isotopic ratio from a Central Mediterranean sediment core. Pre-industrial values are used to train and test two sets of algorithms that are able to forecast the natural variability in δ 18 O over the last 150 yr. These algorithms are based on autoregressive models and neural networks, respectively; they are applied separately to each of the δ 18 O series’ signiﬁcant variability components, rather than to the complete series. The separate components are extracted by singular-spectrum analysis and have narrow-band spectral content, which reduces the forecast error. By comparing the sum of the predicted low-frequency components to its actual values during the Industrial Era, we deduce that the natural contribution to these components of the modern δ


Introduction
Many different archives have been analyzed in the Mediterranean area for the reconstruction of the main physical and chemical parameters characterizing climate over the last millennia. In response to the need for describing natural climate variations and for placing the anthropogenic influence in the right perspective, time series of various lengths and time resolution have been obtained. Just to quote a few examples, we may mention the studies based on vermetid reefs (Silenzi et al., 2004;Sisma-Ventura et al., 2009), on corals (e.g. Montagna et al., 2008), and on marine cores (e.g. Incarbona et al., 2008Incarbona et al., , 2010Piva et al., 2008). However, long series with high resolution (<10 yr) in this region are still rare.
In a previous paper ), we presented a high-resolution record of foraminiferal δ 18 O isotopic ratio that covers the last two millennia. This record was obtained from a shallow-water sediment core drilled in the Central Mediterranean (Gallipoli Terrace in the Gulf of Taranto, Ionian Sea) and was dated with high accuracy by tephroanalysis and radiometric measurements. The δ 18 O series spans the last 2200 yr and was analyzed by singular-spectrum analysis (SSA; see Vautard and Ghil, 1989;Ghil and Vautard, 1991;Plaut et al., 1995;Ghil and Taricco, 1997;Ghil et al., 2002 and references therein). This analysis revealed highly significant decadal, centennial and multicentennial oscillatory components, along with a millennial trend.
The δ 18 O profile shows a steep decrease during the Industrial Era (see gray light curve in Fig. 1). In this paper, we address the problem of evaluating the contribution of natural climate oscillations to this modern variation; pre-industrial δ 18 O variations are used to design and tune algorithms able to forecast the natural variability in the δ 18 O series over the last 150 yr. The comparison between the forecast and the actual δ 18 O signal during the Industrial Era allows one to quantify what percentage of the modern δ 18 O decrease can be attributed to natural vs. anthropogenic causes. Year AD δ 18 O (permil) Fig. 1. Time series of δ 18 O from the GT90/3 core and its components. Gray light curve: raw data; black heavy curve: sum of centennial-to-millennial components + δ 18 O mean value; red heavy curve: trend + δ 18 O mean value. The δ 18 O-axis is reversed, here and in subsequent figures, to agree in tendency with temperature. The SSA results plotted here correspond to a window width of M ≈ 600 yr; they have been shown to be significant at the 99 % level and robust for a wide range of window widths by Taricco et al. (2009 Two independent methods of time series prediction are considered here: (i) autoregressive (AR) models (Box and Jenkins, 1970;Childers, 1978;Montgomery et al., 1990); and (ii) feed-forward neural networks (Hagan et al., 1996;Haykin et al., 1999). Since all prediction methods work better on clean, noise-free signals, the predictions here are not performed directly on the δ 18 O series, but rather on each significant variability component extracted by SSA, as suggested by Penland et al. (1991), Ghil (1992, 1993) and Ghil et al. (2002). The δ 18 O forecast is then obtained by summing the contributions of all components.
The paper is organized as follows. In the next section we present briefly the record and its statistically significant oscillatory components. The prediction methodology is introduced in Sect. 3, while the results and their discussion follow in Sect. 4. Section 5 presents the conclusions.

The δ 18 O record and its spectral analysis
The marine sediment core GT90/3 we use here was extracted from the Gulf of Taranto, Ionian Sea (39 • 45 53 N, 17 • 53 33 E), at a depth of about 200 m, and is 3.57 mlong. Several papers (Bonino et al., 1993;Cini Castagnoli et al., 1990, 1992a,b, 1997, 1998, 1999, 2002a,b, 2005Taricco et al., 2008;Vivaldo et al., 2009) described the dating of the core and previous measurements in the shells of the foraminiferal species Globigerinoides ruber performed on samples from it; these measurements covered the last mil-  lennium. Taricco et al. (2009) extended the measurements of the δ 18 O isotopic ratio down-core to cover the last 2200 yr.
The δ 18 O time series spans 188 BC-1979 AD and it consists of 560 samples with a sampling interval of t = 3.87 yr; the raw data set is shown in Fig. 1 (gray light curve). The SSA of this series revealed highly significant oscillatory components of roughly 600 yr, 350 yr, 200 yr, 125 yr and 11 yr . These oscillations were captured, respectively, by reconstructed components (RCs) 2-3, 4-5, 6-7-8, 9-10 and 11-12; a millennial trend, captured by RC 1, is present as well.
In this paper, we are interested in the low-frequency, centennial-to-millennial variability and therefore we neglect the decadal component (RCs 11-12), as well as the interannual, sample-to-sample variability. The total SSA reconstruction of the low-frequency variability (LFV) of interest here, y LF , is given by where y indicates the mean value of the raw δ 18 O data, and y LF is shown as the heavy black curve in Fig. 1. The trend alone is plotted in Fig. 1 as the heavy red curve, while the individual oscillatory components are shown in Fig. 2. The LFV, as defined in Eq.
(1) and plotted as the heavy black line in Fig. 1, captures roughly 40 % of the total variance of the raw data (light gray curve in the figure). Here, the learning section has length N learn = 180 samples and the test section has a total length of N test = 344 samples; the tests are carried out over a sliding window whose length K is fixed, K = 308 samples, and whose position shifts ahead as the lead increases (red segments). The sliding window provides K pairs of predicted−true values for computing the root-mean-square (rms) forecast errors R(L) at every lead time {L: 1 ≤ L ≤ L max = 36}; see text for details. (a2) Example of iterative one-step-into-the-future forecast producing the first three samples of the sliding window in the position that corresponds to lead L = 10 samples. Input true samples are represented by brown segments. (b) An alternative approach with a longer learning section of N learn = 460 samples, a shorter test section of N test = 64 samples, and no sliding window. In this case, the rms error is evaluated globally over all leads.
δ 18 O values from the early-to-mid 19th century on contain both natural and anthropogenic signals.
Thus, the long pre-1840 section of the series, representing natural δ 18 O variability, with N = 524 samples, is employed to first fit and then test the performance of our two prediction methods, based respectively on (i) AR models, and (ii) neural networks. The prediction is next performed on the 1840-1979 AD interval, comprising L max = 36 samples. As illustrated in Fig. 3, the series is thus divided into three sections: learning, testing, and forecasting. Recall that the sampling interval is t = 3.87 yr ≈ 4 yr and so, for instance, L max t ≈ 140 yr = 1979-1840 yr.

(i) SSA-AR prediction
To be precise, the first method is a combined SSA-AR method, as introduced by Ghil (1992, 1993) and further described in Ghil et al. (2002); it is implemented in the toolkit available at http://www.spectraworks.com/web/ products.html. The basic idea is that the RCs in SSA are data-adaptively filtered signals, each of which is dominated by a narrow spec-tral peak (Penland et al., 1991). Such narrow-band signals allow much more accurate predictions than broad-band ones.
We chose to predict separately each of the five narrowband signals contained in our time series, namely the four oscillatory ones shown in Fig. 2 and the trend shown as the red heavy curve in Fig. 1. The five predictions of RCs 1, 2-3, 4-5, 6-7-8 and 9-10 are then summed to yield the total prediction, y pred .
The SSA-AR method turns out to be most reliable and robust when the order M AR of the autoregressive model fitted to the data is similar to the SSA embedding dimension M SSA , M AR ≈ M SSA . This result follows from the way the SSA lagcovariance matrix and Burg's algorithm AR coefficients are computed in the combined SSA-AR algorithm.
Moreover, it is good practice to use an order M AR that, like the SSA embedding dimension M SSA , is not too large with respect to the length N of the time series, since the variance of the AR-coefficient estimates increases with the order. As recommended by Vautard et al. (1992), we use a section of length N learn = 3M SSA for the learning, with M SSA = 60.
Since the actual prediction is performed with a maximum lead of L max = 36, we first evaluate the prediction 834 S. Alessio et al.: Natural variability and anthropogenic effects in a Central Mediterranean core skill by cross-validating over the test section with the same maximum lead.
Our testing approach is aimed at assessing the forecast skill as a function of lead time, in terms of the root-meansquare (rms) error between predicted and true values, which should be as small as possible, as well as the correlation between forecast and true values, which should be as large as possible. We thus need to verify in the test section several predictions with L steps ahead, up to and including L max . This is achieved by carrying out the iterative scheme shown in Fig. 3a1.
In this approach, the length of the test section is N test = N − 3M SSA = 344 samples. A sliding window of fixed length K = N − 3M SSA − L max = 308 samples shifts ahead, depending on the lead time considered for testing, as illustrated by the red segments in Fig. 3a1. For a value L of the lead, we thus set the initial and final points of the sliding window at x i = 3M SSA + L and x f = 3M SSA + L + K, respectively. For instance, when L = 1, the sliding window that provides the K pairs of predicted-true values for skill assessment starts at the sample x i = 3M SSA + 1 = 181 and ends at the sample x f = 3M SSA + K + 1 = 489 while, when L = L max , the sliding window starts at the sample x i = 3 M SSA + L max = 216 and ends at x f = N = 524, which is the last sample (1840 AD).
For any position of the sliding window, corresponding to a lead time L, each one of the K samples that form the window is predicted using M true values as inputs: samples from x(n + 1 − M) to x(n) are used to forecast x(n + L), by applying an iterative, one-step-into-the-future procedure, according to the equations In system Eq.
(2), we have dropped the subscript AR from M AR for simplicity, and will do so below. When predicting the first sample of the sliding window in its L-th position, only samples out of the learning section are used as inputs, but moving toward the K-th and last window sample, more and more true values from the test section are employed. The rms error R = R(L) thus evaluated quantifies the prediction skill at the particular lead time L. An example of this procedure, for the case L = 10, is shown in Fig. 3a2. The actual forecast for a given RC series is then calculated for the Industrial Era interval 1840-1979 AD by the same iterative procedure. The contributions of all RCs to the LFV in this interval are finally summed to give the total forecast y pred ; the uncertainty in y pred at the various lead times L is computed as the square root of the sum of the squared rms errors of the individual RC forecasts.
In order to check the robustness of the results obtained by the approach described above -in which only a relatively short learning section, with 180 samples, is used -we applied also another AR prediction approach, in which the learning section is longer. As illustrated in Fig. 3b, it includes the first N learn = 460 samples, from ∼200 BC to ∼1600 AD, while the remaining N test = N − N learn = 64 samples, from ∼1600 AD to ∼1840 AD, are used for comparing true values with those predicted, adopting a range of distinct AR orders; no sliding window is used in this alternative approach.
The optimal AR order for the particular component considered is then chosen on the basis of the rms error evaluated globally over all leads. By using this procedure we obtained substantially the same predictions as with the sliding-window approach, which increased our confidence in N learn = 180 samples being sufficient for the AR method to learn our δ 18 O time series' natural LFV.

(ii) Neural-network prediction
For comparison with the AR results, the forecasting of each of the five narrow-band signals contained in our δ 18 O time series was performed also using feed-forward neural networks (Hagan et al., 1996;Haykin et al., 1999). In order to be able to choose the architecture of the network a posteriori, on the basis of algorithm performance, the prediction was tackled by initializing and training a set of feed-forward neural networks that all had a single output layer and a single hidden layer.
We used several numbers I of samples, 2 ≤ I ≤ 12, in the input vector and several numbers H of neurons, 2 ≤ H ≤ 12, in the hidden layer. The ranges of values for I and H were chosen on the basis of preliminary trials performed on all five narrow-band signals. This approach allowed us to choose the network architecture that best reproduces the samples in the test section, which has fixed length N test , as in Fig. 3b.
For this purpose, a MATLAB code was written, employing the functions of its Neural Network Toolbox. The training was performed using the Levenberg-Marquardt algorithm (Rumelhart, 1986;Hagan and Menhaj, 1994). This algorithm is best suited for the type of problem at hand, in which the net is small or of medium size -with less than one hundred weights as an order of magnitude -and the approximation must be very accurate.
The trained network is used to forecast the considered δ 18 O component in the test section, using an iterative, onestep-into-the-future procedure, as done for the SSA-AR prediction. The last I elements of the learning section are used to predict the first sample of the test section. Then, the first predicted sample of the test section is inserted into the next input vector to get the second forecast, and so on, until the 64-steps-into-the-future forecast is produced and the end of the test section is attained; see again (1); SSA-AR prediction in red; and neural-network prediction in blue. See text for details of the prediction methods. Clearly the agreement between the two forecasts and the observed LFV is excellent for a length of time equal to that of the modern, Industrial-Era predictions in the next figure, i.e., out to roughly 1730 AD; it deteriorates only slightly beyond this range, i.e. for the latter part of the test section here.  Fig. 2 and for the sum of RCs: observed LFV in black, cf. Eq. (1); SSA-AR prediction in red; and neural-network prediction in blue. See text for details of the prediction methods. Clearly the agreement between the two forecasts and the observed LFV is excellent for a length of time equal to that of the modern, Industrial-Era predictions in the next figure, i.e. out to roughly 1730 AD; it deteriorates only slightly beyond this range, i.e. for the latter part of the test section here.
The best-performing net is chosen, i.e. the one that has the smallest global rms error over the entire test section. This optimal net is then used for calculating the 36 samples of the post-1840 forecast, which are produced in the same way as the test-section forecasts.
The dependence of the training process of each net on the random choice of the initial guesses for weights and biases suggested running the program a number of times for each narrow-band signal and then averaging the forecasts over all the runs. Finally, the average forecasts for the five narrow-band signals were summed to yield the total average prediction y pred .

(iii) Validation of the methodology
The 36 predicted values {y pred (j ) : j = 1, . . . , 36} of δ 18 O for the interval 1840-1979 AD are thus determined, either (i) by an AR(M) model with order M or (ii) by neural networks. These values represent the extrapolated natural LFV of the δ 18 O record and are then compared with the actual y LF observed after 1840; the latter represents the total LFV, natural plus anthropogenic, for the Industrial Era (see Sect. 4 below). Our comparison of natural and anthropogenic variability thus eliminates higher-frequency contributions, such as year-toyear or, more precisely, sample-to-sample fluctuations.
Before proceeding to the evaluation of the varying contribution of natural and anthropogenic causes to the LFV in δ 18 O during the Industrial Era, we tested our two methods in a pure-prediction mode for a test section that is 64 points (248 yr) long and lies toward the end of the pre-industrial era; to be precise it spans 1592-1840 AD (Fig. 4).
For greater clarity, we only show in Fig. 4 a direct comparison between the two predictions -using the SSA-AR and neural-network approaches -on the one hand, and the LFV of the observed data, on the other, according to the method illustrated in Fig. 3b. The test section in Fig. 4 is almost twice as long as the modern forecast section, while the agreement between the two forecasts, as well as between either of the two and the observed LFV, is clearly excellent for a length of time equal to that of the modern, Industrial-Era predictions; it deteriorates only slightly beyond this range, i.e. for the latter part of the test section.  show the forecast of the five narrow-band signals that compose our δ 18 O record, while Fig. 5f shows their sum y pred . The black heavy curves in all panels represent observed values. In particular, the one in Fig. 5f corresponds to the one shown in Fig. 1, i.e. y LF , but is limited now to the interval ≈1740-1979 AD. Figure 6 compares the total AR forecast over the Industrial Era (1840-1979 AD, red curve), accompanied by its error bounds (pink shading), with the low-pass filtered observed δ 18 O, namely y LF (black curve), put now into the perspective of the entire 2200-yr record. In this figure, we show also the forecast obtained by the neural-network methodology (blue curve). The latter prediction stays within the error bounds of the SSA-AR forecast.

Results and discussion
We notice, in particular, that the natural-variability forecast y pred , whether it is obtained by the AR or by the neuralnetwork approach, reaches δ 18 O values within the Industrial Era, which are as low as those that characterized the Medieval Optimum. The observed y LF values, though, greatly exceed these values, even though our record stops in 1979.
In order to measure the difference between the observed y LF values during the last 150 yr or so and the forecast ones y pred by either prediction method, we define the frac-tion V nat (t) of natural variability by where µ is the mean computed over the pre-industrial section of the time series, i.e. before 1840 AD. The fraction V nat (t) is shown for the AR forecast in Fig. 7 (red curve); the pink shading indicates the range of uncertainty for this fraction, based on the errors in the forecast y pred , using Eq. (3).
In agreement with our assumptions, V nat (t) equals 100 % at the beginning of the Industrial Era (1840 AD in this account), while the anthropogenic contribution remains negligible up to about 1880 AD. Afterwards, the natural portion of the LFV decreases gradually to about 40 % by the end of the 1970s (top of the core). The blue curve represents the same fraction V nat (t) calculated from the neural network forecast; the latter does lie within the uncertainty range of the AR estimates, further confirming our confidence in this range.
? found, for instrumental data over the IE, a natural global climate oscillation of a period 65-70 yr; its exact period, as well as its amplitude, differs depending on the region under consideration. In the case of Eurasia, the period is actually longer than the global one and therefore even more difficult to evaluate with confidence, since it is comparable to the   6. Total δ 18 O SSA-AR forecast y pred (t) (red curve) for the Industrial Era and observed LFV values y LF (t) (black curve) over the entire 2200 yr covered by our record. Both y pred (t) and y LF (t) represent the climatic LFV recorded in the δ 18 O of core GT90/3, namely the ten leading RCs; for consistency, the mean δ 18 O value has been added back in (see Eq. 1). Pink shading: forecast error of the total SSA-AR forecast, based on the five forecasts of the narrowband signals (see Fig. 5); and blue curve: sum of the narrow-bandsignal forecasts obtained by neural-network prediction; see text for details.
total record length of one-and-a-half centuries. It is interesting to note that our 125-yr component (RCs 9-10) forecast over the last 150 yr (bottom panel of our Fig. 2) is in reasonable agreement, within the error bars, with Schlesinger and Ramankutty's (1994) oscillation (see panel 3 in their Fig. 2); in particular, both their oscillation and our forecast share a pronounced maximum (δ 18 O minimum) at 1850 AD and a slightly flatter minimum (δ 18 O maximum) at 1895 AD. Therefore, the natural LFV considered in our study includes in fact, over the IE, most of the variability found by these authors.
Comparing the natural δ 18 O variation y pred and the total one y LF in Fig. 5f allows one to draw several conclusions about temperature variation during the last century. As noted by Taricco et al. (2009), the observed modern variation of about 0.5 ‰ in δ 18 O, if entirely due to temperature effects, would approximately correspond, according to the transfer function of Shackleton and Kennet (1975), to an increase of about 2 • C. This estimate is confirmed by alkenone-derived temperature measurements performed in the same sediments by Versteegh et al. (2007).
According to our prediction, the natural variation in δ 18 O over the last century is of about 0.2 ‰, which corresponds to 0.8 • C. This result suggests that, as early as the end of the 1970s, 60 % of the LFV component of the observed temperature increase was already due to anthropogenic effects.  The high increases in temperature values, observed as well as predicted, indicate a local amplification recorded in the core with respect to the modern variation in Northern Hemisphere temperatures. As argued already by , this amplification is most likely related to the shallow, semi-enclosed nature of the Gulf of Taranto.

Conclusions
The δ 18 O record in the shallow-water sediment core GT90/3 exhibits a marked modern decrease . The remarkable length of this accurately dated δ 18 O series, extending back to 188 BC (Fig. 1), allowed us to separate natural from anthropogenic variability in the modern portion of the record, 1840-1979 AD.
In order to evaluate the contribution of natural variability to the modern decrease in δ 18 O, we fitted autoregressive (AR) models and trained neural networks on the pre-industrial (188 BC-1840 AD) δ 18 O values only. As a noise reduction strategy (Penland et al., 1991), we first applied singular-spectrum analysis (SSA) to extract the statistically significant signals associated with the climate's lowfrequency variability (LFV); these signals include a millennial trend, captured by reconstructed component (RC) 1 (red curve in Fig. 1), and four oscillatory modes, captured by RCs 2-3, 4-5, 6-7-8 and 9-10, respectively (Fig. 2). These five narrow-band signals were then fitted and forecast individually, rather than trying to do so for the original, noisy time series.
Both the AR and neural network algorithms were trained and tested on the pre-1840 part of the record, according to the methods outlined schematically in Fig. 3. The post-1840 forecast values of the narrow-band signals (Fig. 5a-e) were then summed (Fig. 5f) and the sum y pred (t: 1840 ≤ t ≤ 1979) was assumed to represent the natural LFV during the Industrial Era.
The predicted LFV y pred was compared to the total one y LF (Fig. 6), thus allowing us to assess the contribution of natural variability to the LFV of δ 18 O over the last century-and-ahalf. This contribution turned out (Fig. 7) to have decreased steadily from 100 % in the mid-19th century to about 40 % by the end of the 1970s, while the combined high-frequency variability and red-noise background decreased from 65 % in the pre-industrial era to 32 % in the Industrial Era.