2.1 An “accumulation model” of carbonate platform development

Palaeogeographic maps of carbonate platforms spanning the Phanerozoic since the Ordovician in Kiessling et al. (2003) have been used to model carbonate platform accumulation through time. The 31 maps represent the spatial extent of carbonate platforms at the maximum marine transgression within each epoch for the Phanerozoic (Kiessling et al., 2003). Additionally, all instances of carbonate lithologies of Precambrian (4600–541 Ma) and Cambrian age (541–485 Ma) are extracted as point data from the Paleobiology Database (2017) (304 data points; https://paleobiodb.org/, last access: 15 May 2018), the PaleoReefs Database (19 data points; Kiessling and Flugel, 2002) and some well-studied carbonates in northern China (Meng et al., 2011). Considering that the spatial extent of the Precambrian carbonate platforms is poorly constrained, we applied a 200 km buffer to the point data from the databases. All carbonate platform geometries are georeferenced to WGS84 co-ordinates in ArcGIS (ESRI, 2011) to be compatible with GPlates, the open-source plate reconstruction software (Boyden et al., 2011). A GPlates file is created to depict carbonate platform accumulation from 410 Ma to the present day using platforms which have existed since the Precambrian, in order to better capture buried carbonate platforms at the beginning of the model time frame. Creating this model involves stitching together all carbonate platform shapefiles and allowing them to persist until present day. Ages were ascribed based on the Golonka and Kiessling (2002) Phanerozoic timescale. However, owing to new chronostratigraphic data, the geological times ascribed to the existence of carbonate platforms in our model for each epoch are converted to the corresponding times in the latest chronostratigraphic scheme given in the 2016 version of the International Stratigraphic Chart (Cohen et al., 2013). Subsequently, the static carbonate platform shapefiles are reconstructed with the rotations and plate geometries associated with an older plate motion model (Scotese, 2016), similar to the model used by Kiessling et al. (2003). In GPlates, the carbonate platform geometries are assigned plate IDs based upon their overlapping position within the reconstructed continental geometry. Carbonate platforms are rotated to their present-day positions, and plate IDs associated with the Scotese (2016) model were translated into corresponding plate IDs from the Matthews et al. (2016) plate motion model. By obtaining the present-day geometry of the ancient carbonate platforms, attaching these geometries to any other plate motion model becomes straightforward. The evolving carbonate platform model is created by layering each carbonate platform geometry such that it persists to present day. The accumulation model we implement in this study embodies the idea that carbonate platforms accumulate in crustal reservoirs through time.

Figure 1(a) Schematic representation of the ∼ 450 km long cross-profile “whiskers” that are constructed along subduction zones (purple and black), spaced 10 km apart in the direction of subduction. The carbonate polygons (light blue) of the carbonate platform mask have a value of 1 (including those on the continental and oceanic crust), and elsewhere the grid values are zero. All whiskers that intersect, at some point along their length, with the carbonate platforms in the mask will contribute 10 km to the total count of carbonate-intersecting subduction zones. (b) Schematic representations of carbonate-intersecting continental subduction zones. Carbonate platforms become buried over time, forming reservoirs in the crust. Through assimilation and decarbonation reactions, arc magmas interact with upper-crustal carbonates and liberate significantly more CO₂ emissions that at non-intersecting continental arcs.

Download

2.2 Assumptions and limitations of the accumulation model

The accumulation model is an end-member scenario of how carbonate platforms evolve through time. Carbonate accumulation is assumed to accrete on primarily continental lithosphere and persist as mid-crustal carbonate reservoirs from the time of their formation to the present (Fig. 1). Ancient carbonate platforms are known to persist to the present day in surface reservoirs, and can be reconstructed from the geological record as they outcrop, are sampled by drilling or are interpreted from seismic reflection studies (Kiessling et al., 2003). Some crustal carbonate is inevitably eroded or subducted into the deep mantle, however it is improbable that most reservoirs have been destroyed in this way as most carbonate platforms accumulate on continental margins and are not likely to subduct (Lee and Lackey, 2015). Given limited erosion on passive continental margins, except in the cases of major uplift and deformation through mountain-building, it is far more likely that carbonate platform accretion has exceeded their depletion by erosion or consumption through time (Ridgwell and Zeebe, 2005). Moreover, the existence of fossil reef data as far back as the Precambrian (Grotzinger and James, 2000) favours the notion that extensive portions of platforms have been well preserved in continents.

Nevertheless, our model assumes that carbonate platforms have not been significantly depleted by sustained mantle-crust interaction through time. We follow this assumption because there is no way to account for their rate of depletion given the complexity of inter-dependent factors such as platform thickness, heat, pressure, composition and the duration of interaction (Johnston et al., 2011). Accounting for the thickness of platforms and the depletion of reservoirs over time are areas of necessary future developments in the model. We accept that our accumulation model is a simplification of the carbonate reservoir system and may lead to an overestimation of CISZ lengths. However, this may be a reasonable proxy for the substantial amount of CO₂ emitted via diffuse outgassing that is largely unaccounted for in global estimates of carbon emitted at continental arcs (fore-arc, arc and back-arc settings; Keleman and Manning, 2015). Given little evidence of carbon storage in arc crust, these estimates provide important information on, if not an upper-bound estimate of, degassing from CISZs, which are hypothesised to be a significant source of atmospheric CO₂ (Keleman and Manning, 2015). While some Precambrian carbonate platforms have been described (Grotzinger and James, 2000), no effort has yet been made to map their occurrence globally. It is beyond the scope of this work to map the global spatio-temporal patterns of Precambrian carbonate platforms, and thus only those point data captured by the PaleoReefs and Paleobiology databases were incorporated into our accumulation model. While it is expected that accounting for the entire extent of Precambrian platforms would not drastically change the analysis if Precambrian platforms largely lay over the extent of platforms already accounted for, the accumulation model nevertheless represents a degree of underestimation of the Earth's crustal carbonate reservoir.

2.3 Measuring global subduction zone lengths with pyGPlates

We use the open access global plate motion model from Matthews et al. (2016) to analyse the spatio-temporal distribution of subduction zone volcanic arcs in a deep-time tectonic framework and test the hypothesis that subduction zone arc lengths are correlated with atmospheric CO₂ levels. The model spans much of the Phanerozoic as it links the 410–250 Ma Domeier and Torsvik (2014) and the 230–0 Ma Müller et al. (2016) plate motion models. Using open-source and cross-platform plate reconstruction software GPlates (http://gplates.org/download.html, last access: 20 June 2018) and a new pyGPlates Python application programming interface (API; https://www.gplates.org/docs/pygplates/, last access: 20 June 2018), we measure the length of subduction zones at 1 Myr intervals since 410 Ma and track the length of subduction zones that interact with carbonate platforms (CISZ), the length of all subduction zones that do not intersect carbonate platforms (non-CISZ), and total global subduction zone lengths. Using the pyGPlates workflow (Doss et al., 2016), global subduction zone boundaries were extracted from the Matthews et al. (2016) plate motion model, and subduction zone geometries were adjusted to remove any overlapping line segments that would overestimate subduction zone lengths through time (Doss et al., 2016).

2.4 Computing the intersections of carbonate platforms and subduction zones

We test whether subduction zones interacting with carbonate platforms are different from non-intersecting portions of subduction zones in their contribution to global atmospheric CO₂ concentration, and hence we track the lengths of subduction zones that both do and do not intersect with carbonate platforms in the overriding lithosphere. Carbonate platform geometries are converted to polygons using grdmask from the development version of Generic Mapping Tools (GMT; v.5.2.1; Wessel et al., 2015), which create a Boolean-style mask consisting of closed domains with a value of 1 where carbonate platforms existed and 0 elsewhere. We identify polygons proximal to subduction boundaries through time using the exported subduction zone geometries from the pyGPlates workflow. We compute the signed-distance function, positive toward the overriding plate, on the surface of the sphere for all subduction zones and computed where the carbonate platforms lie within +448 km of the trench. The grdtrack tool in GMT is used to create large tracks of 448 km long cross-profiles perpendicular to subduction boundaries at a uniform spacing of 10 km, with “whiskers” pointing in the direction of subduction (Fig. 1a). The whiskers have five evenly spaced nodes to detect intersections with carbonate platform polygons anywhere along the cross-profiles. The whiskers function as a buffer, allowing us to identify the lengths of subduction zone boundaries within ∼ 450 km perpendicular radius of a carbonate platform polygon. Cross-profiles that overlap with carbonate platform polygons demarcate areas where continental subduction zones interact with crustal carbonates (Fig. 1b).

We statistically determine the 448 km long inclusion distance by calculating average arc-trench distances using the present-day “Volcanoes of the World” database maintained by the Smithsonian Institution's Global Volcanism Program (2003). From a sample size of 1023 volcanoes, average arc-trench distances are 287 km with a standard deviation of 161 km, such that one standard deviation greater than the mean arc-trench distance gives 448 km. This captures 84 % of the location of all present volcanic arcs. While flat slab subduction may cause arc-trench distances in excess of 500 km (Espurt et al., 2008), a 448 km long whisker length is likely to capture the majority of carbonate intersections through time.

Data filtering

The modelled results and proxy CO₂ time series were pre-processed for wavelet analysis by detrending and applying a low-pass filter. A moving average filter was applied with a window size of 5 Myr to remove high-frequency signals that may represent noise or other short-term perturbations (Fig. S2). In the time-frequency domain, our focus is on the power in the low-to-mid-frequency spectrum and so the moving average filter was designed to remove periods less than 5 Myr. As a result, the analysis does not capture any geologically rapid (few Myr) changes in CO₂ that may be driven by changes in plate tectonics and carbonate platform-arc interactions.

Rather than comparing estimates from the Accumulation Model with Phanerozoic atmospheric CO₂ models like COPSE (Bergman et al., 2004) or GEOCARBSULF (Berner, 2006), results were compared to the proxy CO₂ record compiled by Foster et al. (2017) because the unresolved uncertainties in model calculations may have rendered comparisons less meaningful than when comparing with proxy records. The multi-proxy CO₂ dataset presented by Foster et al. (2017) is the most up-to-date compilation of CO₂ data at the highest temporal resolution available (0.5 Myr sampling interval), representing their highest probability estimate of CO₂ at each interval and is the proxy CO₂ signal used in this study.

4.1 Continuous wavelet transform (CWT)

The wavelet power spectrum in the CWT shows the spatio-temporal distribution of amplitudes of correlation between the time series and a Morlet wavelet of a specific scale. Areas outside the cone of influence (COI) are not localised in time and are thus unreliable for constraining the duration of periodicity. For this reason, results within the 5 % significance level outside the COI (in white translucent regions of the wavelet power spectra) are ignored.

There is significant wavelet power in the 24–48 Myr waveband (period length) for CISZ lengths, which is intermittent in time yet dominates for time intervals (365–295 Ma, 250–115 Ma, 90–35 Ma; Fig. 2a). Significant wavelet power also exists in the 48–64 Myr waveband of 175–100 Ma and in the 72–128 Myr waveband of 220–120 Ma. In contrast, periodic signals of length ∼ 24–32 and 48–64 Myr have the highest wavelet power and are significant in the non-CISZ time series from ∼ 360–150 Ma (Fig. 2b), and are notably absent from the CISZ time series. Total global subduction zone lengths have two clear and persistent signals through time, primarily in the 24–32 and 48–64 Myr wavebands (Fig. 2c). Significant, high wavelet power in the 16–32 Myr waveband exists in the proxy CO₂ time series of 280–150 Ma and of 125–50 Ma (Fig. 2d). The ∼ 96 Myr waveband has the highest wavelet power at the 5 % significance level in the proxy CO₂ time series of 300–120 Ma (Fig. 2d). There is overlap in periodic behaviour between the CISZ, non-CISZ, global subduction zone and CO₂ proxy data in the 16–32 Myr wavebands.

Figure 2Continuous wavelet transforms (CWT) for four time series: (a) carbonate-intersecting subduction zone (CISZ) lengths, (b) non-CISZ lengths, (c) global subduction zone lengths and (d) the proxy CO₂ data from Foster et al. (2017). Thick black contours designate the 5 % significance level against a red noise background spectrum. The cone of influence (white translucent region) is where edge effects distort the spectrum, meaning that results in this region are not reliable. The colour bar indicates the range of wavelet power in the wavelet power spectrum, with hotter colours corresponding to the maximum peaks in wavelet power. Note the logarithmic scale of the period (Myr) axis.

Download

4.2 Cross-wavelet transform (XWT) and wavelet transform coherence (WTC)

The periodic components common in the subduction zone length time series and the proxy CO₂ record were investigated further using XWT and WTC (Fig. 3). The XWT highlights areas of common high wavelet power between two CWTs, whereas the WTC finds the correlation between two CWTs in time, highlighting persistent and significant periodic signals even if cross-wavelet power is low (Grinsted et al., 2004). Phase-locked arrows within the 95 % confidence contours in the WTC (arrows that do not change direction through time) are the strongest indicator of a persistent periodic signal. Generally, the WTC analysis between subduction zone lengths and the proxy CO₂ record does not highlight significant and prolonged coherence through time, but rather, highlights time intervals where periodic behaviour was significant and stable.

While a prolonged interval of high joint power in the XWT spectra between CISZ lengths and proxy CO₂ exists in the 24–48 Myr waveband (Fig. 3a), this is not found to be coherent through time in the WTC (Fig. 4d). However, the ∼ 96 Myr waveband is coherent in 200–100 Ma and likewise, the 32–48 Myr waveband is coherent in 60–40 Ma (Fig. 3d). Portions of these regions lie outside the COI and are thus cannot be well localised in time. Other areas of high coherence are not linked to significant cross-spectral power in the XWT.

East-pointing arrows that are phase-locked in the WTC in the ∼ 96 Myr waveband indicate in-phase peaks, such that local maxima in CISZ lengths and palaeo-atmospheric CO₂ are aligned over the 200–100 Ma time interval (Fig. 3d). Between 60 and 40 Ma, coherent behaviour between CISZ lengths and the proxy CO₂ record is characterised by a 32–48 Myr periodicity and CISZ lengths lead CO₂ peaks by ∼ 90^∘ (∼ 8–12 Myr; Fig. 3d). Over a similar time interval, coherent behaviour is found between the proxy CO₂ record and global subduction zones. Between ∼ 60 and 20 Ma in the 16–24 Myr waveband, CO₂ peaks lag by 5–8 Myr behind peaks in global subduction zone lengths (Fig. 3f). However, cross-wavelet power is stronger between the CO₂ record and CISZ lengths compared to global subduction zone lengths over the 60–20 Ma interval (Fig. 3a, c).

Non-CISZ lengths and proxy CO₂ record share broad areas of high common power (Fig. 3b), yet the most significant region in the WTC is the 12–16 Myr waveband of 275–225 Ma (Fig. 3e), corresponding to a significant area of the XWT. In this region, north-pointing phase arrows indicate a 90^∘ (∼ 3–4 Myr) lag of peaks in the proxy CO₂ record compared to peaks in non-CISZ lengths. Similarly, peaks in the proxy CO₂ record lag global subduction zone lengths by 3–4 Myr in a cyclic period of 12–16 Myr (Fig. 3e). South-pointing phase-locked arrows indicate that CO₂ peaks are leading peaks in non-CISZ lengths (48–64 Myr waveband, 315–185 Ma) and west-pointing arrows indicate anti-phase behaviour (where peaks in one time series coincide with troughs of another time series; 16–32 Myr, 210–150 Ma), both of which do not support causal relationships.

Figure 3The cross-wavelet transforms (XWT; top) and wavelet transform coherence (WTC; bottom) for (a, d) CISZ, (b, e) non-CISZ and (c, f) global subduction zone lengths respectively and the proxy CO₂ record (Foster et al., 2017). Thick black contours designate the 5 % significance level against a red noise spectrum. The white translucent region designates the cone of influence (COI). The colour bar indicates the magnitude of cross-spectral power for the XWT, and for the WTC it represents the significance level of the Monte Carlo test. The arrows indicate the phase relationship of the two time series in time-frequency space, where east-pointing arrows indicate in-phase behaviour, west-pointing arrows indicate anti-phase behaviour, north-pointing arrows indicate that the peaks in subduction zone lengths leads the CO₂ peaks and south-pointing arrows indicating that CO₂ peaks lead peaks in subduction zone lengths.

Download

The XWT shows that global subduction zone lengths and the proxy CO₂ record share regions of high joint wavelet power (Fig. 3c). Nevertheless, extended horizontal regions in the WTC with south-pointing arrows indicate that CO₂ peaks lead subduction zone lengths (48–64 Myr waveband, 315–185 Ma) and areas of coherence coincide with areas of low cross-spectral power (wavebands < 16 and 16–24 Myr, 380–340 Ma), undermining the potential of causal relationships (Fig. 3f).

Figure 4Graph of total global lengths of subduction zones (black) compared with lengths of all subduction zones that intersect with buried crustal carbonate platforms (magenta, thick) and those subduction zone lengths that do not intersect crustal carbonate platforms (dark magenta, thin). The proxy CO₂ record (blue) from Foster et al. (2017) with upper and lower bounds (blue envelope) is displayed for comparison.

Download

4.3 Comparing carbonate-intersecting arc length and palaeo-CO₂ trends

At present day, our estimates of CISZ and global subduction zone lengths reach approximately 35 300 and 62 000 km respectively. CISZ lengths show an increasing trend over time, especially from 220 Ma to present when lengths reliable reach above 10 000 km globally. From 250 Ma, accumulated carbonate platforms consistently intersect ∼ 70 % of total global subduction zone lengths (Fig. 4). Coherence analysis suggests that CISZ lengths and palaeo-atmospheric CO₂ peaks may not be well correlated prior to 200 Ma. Indeed, CISZ lengths make up a small proportion of global subduction zones prior to 200 Ma, beginning at 8600 km at 410 Ma and increasing to over 14 000 km only after 200 Ma (Fig. 4). During some time windows, periodic behaviour may be correlated. These two areas of significance highlighted by the XWT and WTC analysis are between 200 and 100 Ma, and between 60 and 40 Ma (Fig. 3a, d). Between 200 and 100 Ma, WTC highlights linked periodic behaviour between CISZ lengths and the proxy CO₂ record for a period length of ∼ 96 Myr (Fig. 3d). CISZ lengths increase non-uniformly over this period from ∼ 14 300 to ∼ 23 800 km, whereas palaeo-atmospheric CO₂ levels drop from ∼ 1700 parts per million (ppm) to ∼ 860 ppm. Over the 60–40 Ma period, where WTC analysis found that CISZ lengths lead CO₂ peaks by ∼ 8–12 Myr (Fig. 3d), CISZ lengths show small-scale oscillations between 28 400 and 33 700 km compared with broader, defined peaks in atmospheric CO₂ (Fig. 4). Note that WTC analysis also found coherence with the proxy CO₂ and global subduction zone lengths over a similar period (∼ 60–20 Ma), where peaks in global subduction zone lengths lead peaks in CO₂ by ∼ 5–8 Myr in a 16–24 Myr cycle (Fig. 3f).

At 258 Ma, non-CISZ lengths reach a peak of 86 424 km as global subduction zone lengths reach 95 500 km, followed by palaeo-atmospheric CO₂ reaching a peak of 880 ppm at 252 Ma. This result coheres with cyclic behaviour found in the WTC, where peaks in non-CISZ and global subduction zone lengths precede local CO₂ maxima by 3–4 Myr in the 275–225 Ma time frame (Fig. 3e). In the same interval, a peak at 230 Ma in total subduction zone lengths and non-CISZ lengths (73 700 and 63 600 km respectively) precede a peak at 226 Ma in the proxy CO₂ record (∼ 2100 ppm).

5.1 Present-day subduction zone estimates

At present day, our estimates of CISZ lengths are approximately 2.4 times greater than the recent calculations of present-day global continental arc lengths by Cao et al. (2017; 62 000 and 15 000 km respectively). Cao et al. (2017) focussed on measuring the extent of continental volcanic arcs, whereas our investigation encompasses the total segment lengths of all subduction zones, including those that intersect carbonate platforms. Therefore, our analysis tends to capture a broader range of CO₂ degassing mechanisms along subduction zones, including diffuse outgassing from fore-arcs, back-arcs and volcanic arcs along faults, fractures and fissures, in addition to arc volcanoes. This is because melt is ubiquitous along subduction zones, and there is very little carbon storage along arcs (Keleman and Manning, 2015). For this reason, our approach differs from that of Cao et al. (2017), as we are attempting to capture additional mechanisms of CO₂ escape along continental subduction zones, as opposed to only those directly associated with volcanic arcs, and our estimates may be regarded as an upper bound on total CO₂ emissions at subduction zones.

5.2 Implications for Palaeogene climate

Wavelet analysis also reveals coherent behaviour in the 60–40 Ma interval where local maxima in CISZ lengths consistently precede local maxima in palaeo-atmospheric CO₂ by ∼ 8–12 Myr in a 32–48 Myr cycle. Between 60 and 40 Ma, subduction zone lengths also were found to have local maxima preceding local maxima in palaeo-atmospheric CO₂ by ∼ 5–8 Myr in a 16–24 Myr cycle. However, cross-wavelet power is stronger between CISZ lengths and the CO₂ record, suggesting that CISZ volcanism on a global scale had marginally more effect on climate 60–40 Ma than global subduction zone volcanism. The results may indicate that increases in global subduction zone lengths and CISZ lengths may have contributed to the Early Eocene climatic optimum (EECO) in 53–51 Ma, where temperatures were 14 ± 3 ^∘C warmer than the pre-industrial period (Caballero and Huber, 2013). While climate interactions and feedbacks leading to the greenhouse state remain uncertain, climate sensitivity to CO₂ forcing is likely (Anagnostou et al., 2016).

The termination of arc volcanism particularly along the Sundaland margin, where a large area of carbonate platforms are buried, contributed ∼ 3000 km to the global reduction in CISZ lengths by ∼ 7000 km between 75 and 65 Ma (Fig. 5e, f). The resumption of subduction along the Sunda–Java trench at 65 Ma described by Zahirovic et al. (2016) and the emplacement of a set of carbonate platforms in the lower latitudes at 63 Ma explain the modelled results of an increase in CISZ lengths at 63 Ma which persist until present day (Fig. 5f). The EECO occurred approximately 10 Myr later, supporting the 8–12 Myr lag in peaks of a ∼ 32–48 Myr periodicity found by the CISZ-CO₂ XWT and WTC (Fig. 3a).

Our analysis shows a concomitant increase in global subduction zone lengths (Fig. 4), which cannot be ruled out as an explanatory source of climate forcing. However, modern volcanic CO₂ output from arcs in Indonesia, Papua New Guinea, parts of the Andes and Italian magmatic province are strongly influenced by carbonate assimilation (Mason et al., 2017) and volcanoes in those regions currently contribute significantly to global atmospheric CO₂ flux (Carter and Dasgupta, 2015). Assuming a similar subduction style, this result gives plausibility to the idea that volcanic activity along the Sunda–Java trench during the Palaeogene also significantly contributed to global CO₂ output.

Previous studies have suggested that global arc CO₂ contributed to the baseline warm climate of the Late Cretaceous and late Palaeogene, due to the greater size of carbonate reservoirs in continents and the increased continental arc lengths (Lee et al., 2013; Carter and Dasgupta, 2015; McKenzie et al., 2016; Cao et al., 2017). Our results agree well with the hypothesis, as a relative increase in carbonate-intersecting subduction zone lengths has been linked to a peak in atmospheric CO₂ in 60–40 Ma, contributing to enhanced CO₂ degassing. However, it cannot be determined from this study whether the increase in global subduction zone lengths is more important than the relative increase in the lengths that intersect carbonate platforms.

Figure 5Plate reconstructions with plate boundaries (black), subduction zones (purple) and distributions of carbonate platforms. Carbonate platforms are distributed according to the accumulation model at key intervals from 400 Ma to present day. Colour bar corresponds to the time during which the carbonate platforms were actively developing in the crust. In addition, we set Precambrian and early Phanerozoic carbonate occurrences to appear at the beginning of the model time frame, as the reconstructions do not yet extend beyond Devonian times.

Download

5.3 Limitations of wavelet analysis

Wavelet analysis is inherently sensitive to the shape of signals, and as such, adjustment of filtering and processing techniques introduces variations in the results. A 5 Myr moving average window filter was chosen to remove high-frequency noise, yet it may have moved trends forward in time, which ultimately changes whether peaks overlap or not and adds uncertainty to results. Peaks in the XWT were robust in that they were found to appear in XWT plots regardless of the intensity of smoothing, yet the significance contours differed strongly due to smoothing. For instance, with a longer low-bandpass filter, small-scale oscillations < 5 Myr were removed, reducing areas of significant power in the shorter wavebands. Similarly, when experimenting with a high-bandpass filter, significant regions were directed away from the longer periodicities of > 64 Myr to shorter periodicities. The low-bandpass filter was chosen to remove the noise for smaller periodicities where fluctuations were less likely to represent trends instead possibly represent artefacts of interpolation, and as a result, this study cannot evaluate short-term relationships between climate and subduction zone volcanism.

Our analysis technique cannot quantify the error that propagates from uncertainty in the proxy CO₂ record and subduction zone length signal through to the correlation uncertainty. Firstly, the Foster et al. (2017) proxy CO₂ record is a maximum probability estimate of atmospheric palaeo-CO₂ from a wider 95 % confidence envelope. While the uncertainties appear to be time-localised, the uncertainty in CO₂ concentrations in the proxy record could not be captured by our method of extracting only one time series for analysis. Similarly, the Matthews et al. (2016) tectonic plate model represents one interpretation of geological data which becomes sparse for deeper geological time, especially for pre-Pangaea times (before ∼ 250–200 Ma). As there are no comparable models, we cannot quantify uncertainty in our subduction length signals. Therefore, we acknowledge that our comparison of two time series with large uncertainties are insufficient in constraining periodic behaviour and future research should involve comparing our estimates with other models and data.

XWT analysis describes areas with high power which can be misleading when the XWT spectrum is not normalised. For instance, areas of high joint power were found where relatively flat peaks in the CISZ length data occurred simply because of the large amplitude of peaks in the proxy CO₂ data. The vast differences in magnitude do not necessarily constitute a causative relationship. WTC analysis does not present this problem because the spectrum is normalised, and hence periodicities with strong power were not necessarily found to be coherent (Maraun and Kurths, 2004). The combination of both measures somewhat circumvents the problem of misleading peaks in the XWT, but not entirely. Given that wavelet analysis is very sensitive to magnitude and timing of peaks, the fact that the uncertainty in our modelled estimates of subduction zone lengths are not well constrained introduces variability in correlative relationships. Wavelet analysis, like other signal processing techniques, is limited by the length of the time series data investigated. Periods of greater than half the time series length cannot be examined using this method (Winder, 2002). Therefore, very long oscillations on the scale of 250 Myr or more could not be derived from our data. One reason for examining patterns on 250 Ma periods is due to the very long-term nature of supercontinent cycles. The accretion of supercontinents such as Pangaea and Rodinia is associated with shortened continental subduction zone lengths, whilst continental break-up and dispersal initiates subduction on the edges of continents, and thus increases continental subduction zone lengths (Lenardic et al., 2011; Lee et al., 2016). These changing modes of supercontinent break-up or accretion have the potential to influence climate. Given that supercontinent cycles exceed 200 Myr in length and exert a large effect on the length of continental arcs (Lee et al., 2013), applying this analysis to plate reconstructions that capture multiple supercontinent cycles (e.g. Merdith et al., 2017) to investigate whether such long-term linkages exist between atmospheric CO₂ and global subduction zone emissions.

5.4 Limitations of model assumptions

The extent to which our models are representative of subduction zone CO₂ emissions is limited by some issues with the subduction zone quantification approach. The subduction zone modelling attempted to approximate atmospheric CO₂ flux at subduction zones without directly measuring it by assuming a constant CO₂ flux and a unit correspondence of subduction zone lengths to volcanoes, vents or fissures. Both assumptions are problematic. There are several dynamic factors that govern subduction zone CO₂ flux such as convergence rates (van der Meer et al., 2014), the composition and volume of subducted sediments (Rea and Ruff, 1996), convergence obliquity (Kerrick and Connolly, 2001), slab volume flux (Fischer, 2008), the relative contribution of metamorphic decarbonation in the crust and in the mantle wedge (Keleman and Manning, 2015) and the efficiency of decarbonation (Johnston et al., 2011). The model does not account for these factors and assumes that the parameters of CO₂ flux have remained constant in time. In addition, the number of volcanoes or vents per unit length of subduction zones is neither constant through time nor correlated, at least for the circum-Pacific arc and Central American arc (Kerrick, 2001). Furthermore, the accumulation model has no way to account for carbonate platform thicknesses, and thus cannot account for the realistic depletion of crustal carbonate reservoirs through time. The inability of the model to incorporate the complexity of depleting carbonate platforms may lead to sampling bias during continental collisions. This highlights an area of future improvement, including the utilisation of subduction zone geodynamic–geochemical modelling (e.g. Gonzales et al., 2016) by tracking different CO₂ reservoirs in a geodynamic context.

Furthermore, uncertainty cannot currently be measured using our GPlates modelling approach. Uncertainty propagates through the plate motion model over deeper reaches of geological time, such that it is currently impossible to quantify uncertainty of all possible configurations of subduction zones, especially in the pre-Pangaea time frame. Moreover, there is much greater uncertainty in the dynamic mechanisms of carbon cycling along subduction zones, and quantification of all geochemical decarbonation reactions is still in infancy (Gonzales et al., 2016). For these reasons, formal uncertainties cannot yet be quantified for models of subduction zone carbon cycling. Taking all limitations into consideration, the accumulation model still provides a first-order approximation of changing subduction zone lengths through time and the interaction with (buried) carbonate platforms as an indicator of CO₂ degassing, especially during the Palaeogene.