Integrating palaeoclimate time series with rich metadata for uncertainty modelling: strategy and documentation of the PALMOD 130k marine palaeoclimate data synthesis

Palaeoclimate data hold the unique promise of providing a long-term perspective on climate change and 15 as such can serve as an important benchmark for climate models. However, palaeoclimate data have generally been archived with insufficient standardisation and metadata to allow for transparent and consistent uncertainty assessment in an automated way. Thanks to improved computation capacity, transient palaeoclimate simulations are now possible, calling for data products containing multiparameter time series rather than information on a single parameter for a single time slice. Efforts are 20 underway to simulate a complete glacial-interglacial cycle using general circulation models (palmod.de) and to confront these simulations with palaeoclimate data, we have compiled a multi-parameter marine palaeoclimate data synthesis that contains time series spanning 0 to 130,000 years ago. We present the first version of the data product that focuses exclusively on time series for which a robust chronology based on benthic foraminifera δ18O and radiocarbon dating is available. The product contains 896 time 25 series of eight palaeoclimate parameters from 143 individual sites, each associated with rich metadata, age-depth model ensembles and information to refine and update the chronologies. This version contains 205 time series of benthic foraminifera δ18O, 169 of benthic foraminifera δ13C, 131 of seawater temperature, 174 and 119 of planktonic foraminifera δ18O and δ13C and 44, 38 and 16 of carbonate, organic carbon and biogenic silica content, respectively. The data product is available in three formats (R, 30 LiPD and netCDF) facilitating use across different software and operating systems and can be downloaded at https://doi.pangaea.de/10.1594/PANGAEA.908831 (Jonkers et al., 2019). This data descriptor presents our data synthesis strategy and describes the contents and format of the data product in detail. It ends with a set of recommendations for data archiving.


35
Global climate has varied dramatically over the last glacial-interglacial cycle. Since the previous interglacial (approximately 130,000 years ago) the earth has slowly cooled until the Last Glacial Maximum (LGM; approximately 21,000 years ago). This cooling was associated with the growth of massive ice sheets in North America and Eurasia, leading to a sea level drop of about 120 metres (Waelbroeck et al., 2002) and pronounced climate variability on millennial time scales (Voelker and et al., 2002). Since the 40 LGM, the Earth has warmed rapidly until the onset of the current relatively stable warm period, the Holocene (Shakun et al., 2012). The ultimate cause of the large-scale variations in the Earth's climate are changes in the orbit of the Earth around the Sun (Hays et al., 1976). However, complex feedback and non-linear mechanisms, involving ocean (atmosphere, cryosphere) circulation and biogeochemical cycles, are required to explain how slow changes in the orbital configuration led to the observed evolution of 45 global climate, and how these processes led to the manifestation of abrupt climate change.
For these reasons the last glacial-interglacial cycle has been a key target for palaeoclimate modelling.
Initially this only involved equilibrium simulations for key time slices, such as the LGM, or transient simulations for short periods, such as the last millennium. The motivation to simulate past climate states 50 is given by the possibility for palaeoclimate data to serve as a benchmark for the models. Indeed, this possibility contributed to the development of large palaeodata syntheses (CLIMAP project members, 1981;MARGO project, 2009). The time-slice modelling approach is still being pursued, for example in phase 4 of the Paleoclimate Intercomparison Project (PMIP), four of the five target intervals fall into the time frame of the last glacial cycle (Kageyama et al., 2018). However, with increasing computing power, 55 the focus is now shifting towards transient climate simulations (Liu et al., 2009;Latif et al., 2016) and the simulation of the last deglaciation is now also considered in the PMIP protocol (Ivanovic et al., 2016). This development calls for a different type of palaeodata synthesis, with focus on time series rather than time slices. Time series of climate data are needed to evaluate aspects of transient simulations that are 60 not available in equilibrium simulations, such as rates of change, phase relationships and spectral properties of climate variability. It is also clear that an evaluation in multi-parameter space using different aspects of the climate system and multiple proxies will be more powerful and diagnostic (Kurahashi-Nakamura et al., 2017), calling for multi-parameter synthesis products.
3 chronological uncertainty. The latter is particularly relevant for the comparison of transient climate change and chronological uncertainty thus requires a comprehensive treatment in data syntheses of palaeoclimate time series.
Accounting for proxy uncertainties in a comprehensive and transparent manner requires not only expert 75 knowledge, but also the availability of extensive metadata in addition to the proxy data. However, due to a lack of standardisation and inconsistent archiving of metadata, synthesizing palaeoclimate data in a way that allows robust uncertainty assessment remains challenging and time consuming. Efforts are underway to alleviate these challenges. The largest palaeoclimate data repositories (World Data Service for Paleoclimatology, operated by the national centers for environmental information (NCEI) at NOAA 80 and PANGAEA) are both striving for more standardisation and to store data in (more) machine-readable format. In addition, standardisation is progressing through the use of existing data formats from other communities (netCDF; Langner and Mulitza, 2019) as well as the implementation of new data formats specifically targeted to palaeodata (Linked PaleoData (LiPD); McKay and Emile-Geay, 2016). At the same time there is ongoing discussion on data and metadata requirements and standards (Khider et al., 2019). 85 Traceability of datasets is also improved through data citations, ensuring not only that data producers receive proper credit for their work, but also allowing for better linking of different datasets.
Nevertheless, these initiatives are only recently emerging and the majority of the palaeoclimate data remains inconsistently formatted, non-standardised and scattered over various data repositories. The need for synthesis products and documentation of potential synthesis approaches is therefore as large as 90 ever.
Here we present the first version of a new multi-proxy marine palaeoclimate data synthesis that covers the past 130,000 years developed within the German climate modelling initiative PALMOD (Latif et al., 2016). We focus on the ocean as it is a large reservoir of heat and CO2 and allows global coverage with 95 consistent chronological control. This synthesis goes beyond the time frame of many existing multiproxy/parameter data syntheses (PAGES2k Consortium et al., 2017;Routson et al., 2019), expands existing data products that provide long palaeoclimate time series to multiple parameters (Shakun et al., 2012;Marcott et al., 2013;Peterson and Lisiecki, 2018;Snyder, 2016) and is based on a strategy of semiautomated data harvesting (Cartapanis et al., 2016). This version of the synthesis contains data on nine 100 climate-sensitive parameters: benthic and planktonic foraminifera stable oxygen and carbon isotopes, seawater temperature, radiocarbon and bulk sediment carbonate, organic carbon and biogenic silica content.
In this manuscript we describe our synthesis approach, the contents and structure of version 1.0 of the 105 data product, suggestions for its use, plans for future updates and recommendations for archiving new 4 data and retrieving dark data in a way allowing optimal future re-use. The data product is intended to be used to investigate spatio-temporal changes in a multi-parameter domain. Thanks to rich metadata that allow to rigorously quantify reconstruction uncertainties, we also envision that this data product will provide the building blocks for intelligent palaeoclimate data model comparison , for 110 instance through proxy system modelling (Dolman and Laepple, 2018) or data assimilation (Breitkreuz et al., 2019).
The structure of this data descriptor is as follows. Section 2 describes the synthesis strategy, including the data discovery approach, standardisation and age modelling. In section 3 we provide general information 115 on palaeoclimate proxies from marine sediment archives that is used to guide the metadata selection.
Section 4 details the structure of the database and the contents of version 1.0 are outlined in section 5.
The formats of the data product and where it can be accessed are described in section 6 and we discuss future plans, versioning and intended use in section 7. In the last section, 8, we reflect on the data synthesis effort and provide recommendations for data archiving and data rescue. 120

Data synthesis strategy
Our data product focusses on time series from marine sediment archives. A single marine sediment archive (sediment core) can be used for measurements of different parameters, each providing information on different aspects of the environmental conditions at the time of deposition. However, for the purpose of analysis, the various proxy time series must refer to a single age-depth model for the 125 sediment core they are derived from. For this reason, the basis of our synthesis is formed by a collection of sediment cores, each associated with its own age-depth model.
Marine sediments are dated using absolute age controls, where specific layers are dated using, for instance, radiocarbon, tephra or palaeomagnetic properties, and/or relative age controls, where time 130 series are aligned based on hypothesised synchronicity of the changes recorded by some properties of the sediment. A well-established hypothesis-based age modelling approach with a solid theoretical basis is the alignment of benthic foraminifera stable oxygen isotope ratio (δ 18 O) time series (Lisiecki and Raymo, 2005). We thus base our chronological framework on a combination of radiocarbon dates and benthic foraminifera δ 18 O and have selected time series where both parameters are available as the 135 foundation of this data product. This approach of blending absolute and relative age controls is required to provide age-depth models for sediment cores that extend beyond the radiocarbon dating range (~40,000 years). If available, further proxy time series were then added, thus ensuring a common chronology among all proxy time series measured on the same sediment core.

5
We selected palaeoclimate parameters to synthesise following discussion with climate modellers within the PALMOD project. The high priority selection includes both physically and biogeochemically relevant parameters, of which some are based on measurements that can be compared with climate model output using (forward proxy) models (e.g. benthic δ 18 O) and others represent inferred parameters that can be compared with model output more directly, but for which proxy models are still in their infancy 145 (e.g. temperature based on foraminifera Mg/Ca). Parameter selection also considered the expected spatial and temporal coverage of data availability as well as the existence of previous data products. The high priority parameters for which data are presented here are listed in Table 1. If available, raw data were synthesised and in cases where raw data were not available and it was possible to derive the raw data from the inferred palaeoclimate data, raw data were back-calculated. Raw data time series obtained 150 in this way are flagged with a note describing the calculation.
We note that our approach of first building the stratigraphic framework based on radiocarbon dates and benthic foraminifera δ 18 O means that the synthesis is not necessarily comprehensive as it does not include time series where one of the parameters of interest has been measured, but for which the 155 components of the stratigraphic framework are not available. However, at this stage, we opted to include only sediment cores where an age modelling strategy that is consistent and comparable across the entire data product could be achieved.

160
In principle, data synthesis can proceed by expansion or reduction (Fig. 1). The first, more traditional, approach relies on expert knowledge of what data is available and/or on systematic literature search. In this approach the synthesis grows by including more data until sufficient data that meet inclusion criteria are compiled. In this way, a lot of time is spent on discovering and retrieving datasets, and it is possible that valuable, but less exposed data is missed. On the other hand, this approach has a chance of 165 uncovering dark data that is not publicly available (Fig. 1).
The second approach starts from a large and "crude" synthesis of data from public sources and proceeds by weeding out data that do not meet criteria for inclusion in the data product. This approach faces different challenges: making sure that the initial bulk database is comprehensive (efficient data mining) 170 and assuring that the data filtering is efficient (fast and accurate). In contrast to the expansion approach, this reduction approach cannot discover dark data. However, it is more objective (less reliant on expert knowledge), can be automated more easily and focuses on data that is already in the public domain so that no time is lost on finding data that ultimately proves unavailable. This second approach also rewards and encourages good data stewardship. 175 6 In theory, both approaches can lead to a similarly sized/exhaustive synthesis, but they differ in the allocation of effort (Fig. 1). In practice both approaches are often combined, especially towards the end of a synthesis project, when the data product is benchmarked against existing syntheses.

Initial synthesis
We followed the reduction approach and used a semi-automated pipeline to compile data from public sources. Keywords (supplementary information) were used to make lists of URLs of potentially relevant data on pangaea.de and the linked files were then downloaded in bulk (n = 108,239). A slightly different 185 approach was followed for the NCEI archive. Here, all files that were machine readable at the time of download (September 2016, n = 1,925) were obtained from the ftp server (ftp://ftp.ncdc.noaa.gov/pub/data/paleo/paleocean//sediment_files/complete). Custom scripts in R were used to put all data in a common format and merge time series that could be unambiguously assigned to the same core (based on name and x, y, z position). This resulted in a mixture of records that were 190 merged to the same core and those that could not be either because there was only one data file for the core or because of ambiguous labelling. We refer to the locations of these records as "sites". In order to facilitate the analysis (filtering) of the sites, a uniform attribution of the various parameter names had to be developed. Because no standardised names exist for palaeoclimate parameters, the uniform attribution required the development of attribution libraries for each desired palaeoclimate parameter. 195 The initial synthesis contained time series from 38,511 sites.

Data reduction and standardisation of ontologies
The initial synthesis was reduced by removing non-marine sites (using elevation flag) and further constrained by only considering sites where at least one data point of any of the parameters measured in 200 that core fell within the target time frame (disambiguating age units in the synonym library of the category "age") and the site had benthic oxygen isotope data. This resulted in 781 sites. At this stage, no criteria for length or resolution were applied but we prioritised processing time series that we estimated to contain at least 50 data points within the 130,000-years timeframe. Further data processing started with dereplication of the selected sites. This was necessary because no standards exist for the naming of 205 cores and the repositories store data with different renditions of the same core name, sometimes even associated with erroneous geographic coordinates. This process was done manually and proceeded by constructing a list of disambiguated sites through sequential one-by-one comparison. Where a strict synonym was found (different labels for the same core but the same data), only unique data were retained. At this stage, disambiguation of site names was only done for sites that had at least benthic 210 δ 18 O, so time series of other parameters, but associated with inconsistent core labels could have been 7 missed in the synthesis. However, those sites are contained in the initial bulk synthesis and are hence not lost, but will be salvaged in updates of the data product (see section 7).
Further steps required a manual standardisation of the names of the parameters and their attributes 215 (such as species name that was analysed for oxygen isotopes). This was accomplished by deciding on a final, uniform list of parameters and associated metadata and their possible values. Original parameter names were preserved to allow for cross checking. By metadata, we refer to aspects of the individual parameters that were deemed essential to facilitate a meaningful analysis in a palaeoclimatic context, considering potential sources of uncertainty, such as species of foraminifera analysed or calibration 220 equation used for palaeotemperature estimates. The list of metadata is provided in Table 2. The standardisation was accompanied by further dereplication of individual time series that were already associated with the same site name but archived more than once.

Metadata and chronology
225 Subsequently, as far as possible, metadata values were added manually when missing, often by scraping the information from the original publication. Next, all time series from a single site (core) were put on a common depth scale to allow for age modelling. Data that could not be put on a depth scale were excluded from the synthesis. This was the case where parameter values were recorded only against age and where no other data file was available that allowed unambiguous re-assignment to depth. Ambiguity 230 also resulted from the use of multiple (composite) depth scales for the same archive. Finally, chronological data (all absolute markers, including radiocarbon dates and associated metadata) were manually added, where necessary also by consulting the original publications. Throughout the process publication information (digital object identifier (DOI), or if not available full bibliographic details) and data source (URL and/or DOI) were preserved in order to trace the source of the data. This applies both 235 to the source of the individual data files from repositories and to the sources of the metadata and chronological data.

Age modelling
Whereas the initial steps of data discovery and synthesis could rely on published chronology, analysis of 240 the complete dataset will require the development of a common chronological framework. This framework must be constructed in a way that not only allows a consistent method of assignment of ages to depths within each core, but also allows consistent and quantitative assessment of age uncertainty. To this end, we follow an approach that combines absolute ages (radiocarbon ages, tephra layers and palaeomagnetic events) with δ 18 O stratigraphy. As a result, our age models may differ from those 245 reported in the original publication(s). This is not to state that the updated age models are better 8 (constrained), but they are constructed in a way that allows applicability and consistency across the synthesis. The consistent approach allows an assessment of age uncertainty jointly for all records by a Bayesian approach, generating ensembles of sedimentation histories consistent with the available age control points for each core, allowing uncertainty estimates at each depth by considering the distribution 250 of ages given by the ensemble.
With respect to the reporting of the chronology, we follow a transparent approach, preserving the initial age model and providing the new age models as well as all information needed to revise or update the new age models. In the final step, the age information from absolute ages and δ 18 O tie-points is 255 combined and the age model and its uncertainty was assessed in a Bayesian framework using BACON (Blaauw and Christen, 2011). The entire age modelling routine was carried out in PaleoDataView (PDV; Langner and Mulitza, 2019).
To ensure a common chronological framework for all time series in the synthesis, radiocarbon ages were 260 re-calibrated using the IntCal13 curve (Reimer et al., 2013). Since reservoir ages vary in space as well as in time, we used reservoir age estimates based on a comprehensive ocean general circulation model (Butzin et al., 2017) to account for this variability in a physically plausible way. To derive the reservoir age and uncertainty for a measured radiocarbon age PDV (i) extracts all modelled radiocarbon ages from the nearest grid cell in the modelled data set, (ii) finds all modelled radiocarbon ages that are possible within 265 the error of the measured radiocarbon age and (iii) takes the mean and the standard deviation of all corresponding reservoir ages as correction for the measured radiocarbon age. By definition, this approach cannot account for processes affecting the reservoir ages on subgrid spatial scales. Given the relatively coarse resolution of the model, this means that processes such as upwelling are not fully accounted for. In addition, no modelled reservoir age data are available for the Mediterranean and Red 270 Seas, we use the reservoir ages reported by the authors of the original publication and an assumed uncertainty of 100 years for these basins (5 sites). Absolute ages based on North Atlantic tephra layers and palaeomagnetic events were updated and harmonised using Svensson et al. (2008).
In addition, and beyond the 14 C dating realm (~40,000 years) the age models rely on manual tuning of the 275 benthic foraminifera δ 18 O time series from each core to regional benthic foraminifera δ 18 O stacks (Lisiecki and Stern, 2016). Stable isotope stratigraphy in theory provides a range of events to correlate, however, in order not to inflate confidence in the tuned age models and ensure comparability between different cores, in our approach, the tuning was done as far as possible only by matching the position of marine isotope stage boundaries. We updated the age-depth models only for the 0-130,000 years time frame of 280 this synthesis, but data and original age models extending beyond 130,000 are preserved in the data product. To obtain uncertainty on the age control points obtained by δ 18 O tuning, we used the 9 chronological uncertainty of the δ 18 O stacks, as reported by Lisiecki and Stern (2016). Additional uncertainty associated with the identification of the control points in the individual records or with the assumption on synchronicity was ignored, as these are difficult, if not impossible, to quantify. 285 3. Notes on palaeoclimate proxies in marine sediment archives and metadata This synthesis contains climate sensitive proxy data based on measurements on various biological sensors. It is not the intention here to provide a full overview of marine palaeoclimate proxies and their uncertainties (for this, see e.g. Hillaire-Marcel and De Vernal (2007); Moffa-Sanchez et al. (2019)), but the 290 fact that the proxies are based on biological sensors means that they are affected by different ecological bias in addition to observational noise. Basic knowledge of the recording system is therefore essential for the interpretation of the data and may aid to explain differences between proxies for the same climate parameter. These considerations were also essential to choose the range of metadata to be recorded alongside each palaeoclimatic parameter, to allow a proxy-specific assessment of uncertainty. 295 Foraminifera are among the most widely used proxy sensors in palaeoceanography. They are unicellular marine zooplankton. The species used here all build a calcite skeleton that is preserved in the sediment.
Foraminifera can be divided into two main groups: benthic foraminifera living at the seafloor or at shallow depth in the sediment and planktonic foraminifera living in the upper hundreds of metres of the 300 ocean. In the data product, proxies measured on these two groups are clearly distinguished by a benthic or planktonic prefix. The chemical composition of foraminifera reflects environmental conditions of the seawater the organisms calcified in. For the purpose of this data product the parameters of interest are stable oxygen and carbon isotope and Mg/Ca ratios. Stable oxygen isotope ratios in foraminifera calcite reflect a combination of temperature and δ 18 O of seawater (Urey, 1948), which is in turn related to ice 305 volume and salinity. Species-specific calibrations exist to quantitatively link δ 18 Oforaminifera and δ 18 Oseawater to temperature (e.g. Marchitto et al., 2014;Bemis et al., 1998). Stable carbon isotope ratios (δ 13 C) reflect the δ 13 C of the dissolved inorganic carbon in seawater. In particular the benthic foraminifera species Cibicidoides wuellerstorfii generally incorporates δ 13 CDIC without a biological offset and can serve as a tracer of bottom-water δ 13 CDIC, which is commonly used as non-passive circulation tracer (Curry and 310 Oppo, 2005). The δ 13 C of other benthic foraminifera species is generally not indicative of bottom water δ 13 CDIC and the δ 13 C of planktonic foraminifera is also influenced by temperature and carbonate ion concentration, rendering interpretation complicated (Spero et al., 1997).
The Mg/Ca ratio in foraminifera calcite can be used to infer calcification temperature and in combination 315 with δ 18 Oforaminifera, the δ 18 Oseawater (Elderfield and Ganssen, 2000). Similar to stable oxygen isotopes, species-specific calibrations exist to quantitatively reconstruct past temperature from Mg/Ca ratios (Anand et al., 2003;Lear et al., 2002) and whenever when indicated in the original publication, the calibration is included in the metadata. Carbonate system parameters and salinity have a secondary influence on Mg/Ca ratios in foraminifera calcite (Gray et al., 2018). Whereas benthic foraminifera live in 320 a generally stable environment, the near sea surface habitat of planktonic foraminifera shows large seasonal and vertical gradients. Species-specific seasonal and/or depth habitat preferences may therefore leave a considerable imprint on the proxy signal contained in their shells (Jonkers and Kučera, 2017;Mix, 1987). For all proxies based on foraminifera, it is relevant to record the species as well as the number of individuals that were pooled for geochemical analysis. The latter is because the short life span 325 and variable habitat of foraminifera species causes large variability among individuals. Planktonic foraminifera shell size may for several reason also affect their chemistry (Jonkers et al., 2013;Friedrich et al., 2012) as well as their assemblage composition (Al-Sabouni et al., 2007). Therefore, the size fraction of the analysed shells was included in the metadata whenever this information was available.

330
Besides planktonic foraminifera Mg/Ca ratios, the U K' 37 unsaturation index can provide information about near sea surface temperature. The U K' 37 index is based on the relative degree of unsaturation of C37 alkenones, which is linearly related to temperature (Prahl et al., 1988). Alkenones are produced by coccolithophores, marine phytoplankton living in the photic zone. The production of alkenones is in many regions not constant during the year, thus potentially causing a seasonal recording bias in the U K' 37 335 temperature proxy (Rosell-Melé and Prahl, 2013). Several calibrations exist that relate the index to sea surface temperature and if the calibration was mentioned in the original publication it was preserved in the metadata.
A large proportion of the temperature estimates in this data product are based on microfossil (planktonic 340 foraminifera, diatoms, radiolaria, dinoflagellate cysts) assemblages. These reconstructions are based on a statistical relationship between species assemblages and temperature (Imbrie and Kipp, 1971). In theory, microfossil assemblages can be used to reconstruct temperatures of different seasons, or different environmental parameters from the same assemblage. However, it is not always clear that such reconstructions are truly independent (Telford and Birks, 2011). Several different methods exist to relate 345 fossil assemblages to temperature and researchers often apply more than a single method in their reconstructions to increase confidence (Kucera et al., 2005). When available, these different reconstructions are included in the data product.
The bulk sediment data (CaCO3, TOC and BSi) form a category of their own. They are not proxies in the 350 strict sense, but properties of the sediment that reflect a combination of export productivity, sedimentation and preservation. However, they can provide crucial information about the ocean-climate system, in particular on biogeochemical cycles (Cartapanis et al., 2016). With the advent of explicit sediment modules in climate models (Heinze et al., 1999;Kurahashi-Nakamura et al., 2019), sediment composition can also be directly compared with model output and potentially provide additional 355 constraints on the simulations.

Structure of the database
Following the data synthesis strategy outlined above, we generated a first data product for time series of eight parameters in sediment cores with radiocarbon and benthic δ 18 O stratigraphy. Following the logic of our approach, the synthesis is organised by the physical object from which the records were extracted 360 (cores), here called site, to account for inclusion of records from spliced cores.
Each site in the data product has information on seven different themes (Fig. 2): 1. Geographic data contains the site name, latitude and longitude (in decimals N and E) and elevation/water depth (in m) and possible notes that are relevant to the site/core as a whole. All fields except notes are essential and always included. 365 2. Metadata includes the original parameter name as given in the online data file, a standardised parameter name, parameter type (measured or inferred), unit, an estimate of analytical error as determined from repeat measurements of a standard and of reproducibility as determined by repeat measurements on samples. For measured parameters information on the instrument and laboratory is given. All metadata terms are listed in Table 2 and an overview of the standardised 370 parameter names is provided in Table 3. 3. Chronology data contains raw data on absolute age control points used for age modelling. This includes depth, radiocarbon ages including their uncertainty, dated material, laboratory codes, but also calendar ages of tephra layers and palaeomagnetic events. Age control points not used in the age model by the authors of the original publication(s) are indicated and if available the 375 source (DOI or URL) of the data is shown in addition to the original publication DOI. A complete list of chronology data terms is given in Table 4. 4. The actual time series data are provided on a common depth scale. Original age models are preserved alongside the data time series as they may differ for different time series from the same site. The data also contains information on sample number/label (mainly for 380 DSDP/ODP/IODP cores) for spliced records and sample specific notes.
5. The revised age model contains ages for depths bracketed by age control points (absolute and relative). Mean and median ages are given as well as an uncertainty range (2.5 and 97.5 percentiles) based on the full suite of age model ensembles. No attempts were made to extrapolate the age models beyond the tie points in the 130,000-year time frame, so original age 385 models may extend to either side. 12 6. The BACON data contains all information to reproduce the revised age model. Besides the 14 C and absolute age control points, this includes the tie points for the alignment, the alignment target and all parameters used to construct the age-depth model. 7. Age ensembles are provided for further assessment of chronological uncertainty. In order to 390 keep file sizes manageable 1,000 randomly selected age model ensembles are preserved.

PALMOD 130k marine palaeoclimate data synthesis v1.0 contents
The data product contains 896 time series of the palaeoclimate parameters listed in Table 1 from 143 sites ( Table 5). By design all sites have both benthic stable oxygen isotope and radiocarbon data. The 395 majority of the sites are close to the continents and in the northern hemisphere, with a concentration in the North Atlantic Ocean (Fig. 3). This reflects both research attention as well as the challenges of obtaining sediment cores with high accumulation rates and well-preserved foraminifera. The effect of research focus is also visible in the temporal coverage of the time series, where every parameter is characterised by a clear maximum around the last glacial termination ca. 15 ka BP (Fig. 4). The median 400 resolution of the time series varies by two orders of magnitude, but is generally better than 1 sample per 1,000 years and fairly similar among the different parameters (Fig. 5). The updated age models are based on chronological control points (tie points) from radiocarbon dating, absolutely dated layers (using tephra and/or palaeomagnetic event stratigraphy) as well as alignment to the regional benthic δ 18 O stacks. The majority of the time series have a chronological control point at least every 5,000 years (Fig. 6). Taken 405 together, the coverage in space, time and across parameters indicate that the PALMOD marine palaeoclimate data product allows for analysis of palaeoclimate on a supra-regional scale over the entire 130,000-year time frame.
This data product builds upon previous syntheses. Virtually all of the sites are also part of the benthic 410 foraminifera δ 18 O and δ 13 C compilations of Lisiecki and Stern (2016) and Peterson and Lisiecki (2018). Our synthesis however, also includes data on other palaeoclimate parameters and contains more metadata and information on the age-depth models. Some of the planktonic foraminifera δ 18 O and Mg/Ca time series in the PALMOD 130k data product are also included in the iso2k synthesis effort (Konecky et al., 2018) and a number of sea surface temperature time series are also part of the forthcoming 415 temperature12k synthesis (Kaufman et al., submitted).

Data formats and access
We provide the data products in three different formats in order to facilitate access and analysis using different software and across operating systems. Given the structure of the formats, each representation 13 is slightly different in its level of metadata detail and the way metadata and data are stored. The 420 differences are described below.
• Since the data product was built using R the data product is presented in R-readable RDS files that contain for each site a list with data for each theme (section 4). This is the format that is most complete, yet in the interest of memory space it preserves a random selection of 1,000 age models from the larger ensemble produced using BACON. In this format, all data and metadata 425 for each site are contained within a single file. Sample scripts (https://github.com/lukasjonkers/PALMODutils) allow the user to extract a quick overview of the contents of the data product similar to Table 5, but with additional information on the temporal range, resolution and age control of the time series. Additional code is available to query the data product by parameter, parameter detail, sensor species, temporal range, resolution and age 430 control.
• The data product is also provided in the LiPD format, which is built around JSON-LD and csv formats and is widely readable across different platforms. As for RDS, all data for each site is presented in a single file. Utilities to interact with LiPD files in R, Matlab and Python are available at https://github.com/nickmckay/LiPD-utilities. 435 • Finally, the data is also provided in netCDF format in a way that allows reading in PDV. This means that a single site has separate files for each individual palaeoclimate parameter as well as for the age model. In this format, some metadata is stored as concatenated strings, rather than easily searchable attributes. The netCDF format allows however to store the full suite of age model ensembles without excessive file sizes. The PDV software to read and process the data can 440 be downloaded at https://www.marum.de/en/Stefan-Mulitza/PaleoDataView.html.
The data product is available for download at https://doi.pangaea.de/10.1594/PANGAEA.908831 . We encourage users of the data product also to cite the primary source of the data when using (individual time series of) this product.

445
To increase the spatio-temporal coverage over the entire 130,000-year time frame of the database, updates of this data product will first aim for quantitative growth of the database by adding more time series with chronological control based on benthic foraminifera δ 18 O and absolute age control points other than 14 C. If available these updates will also include the parameters listed in Table 1. They will be named using the counter following the first decimal separator. The structure of the data product is 450 designed to be flexible, allowing for the addition of different metadata fields and parameters. Further updates that include new parameters and/or require a new age modelling approach (i.e. no benthic δ 18 O alignment) will be named using the counter before the first decimal separator. Any updates to add or correct (meta)data to an existing version and that do not increase the number of sites will be indicated using a counter following the second decimal separator. 455

Lesson learned: recommendations for data archiving
Data re-use and sharing are both made easier when data are archived in a standardised manner. Even though a large amount of palaeoceanographic data is publicly available, metadata to facilitate interpretation of the raw or inferred palaeodata is often not, or only partially, made available and needs to be obtained from the original publication. Synthesis efforts therefore still require a lot of time and 460 effort to find, compile and standardise data and metadata. Only recently has the palaeoclimate community started to discuss data archiving standards (Khider et al., 2019). However, implementation of the proposed Paleoclimate Community reporTing Standard (PaCTS 1.0) will only affect new uploads to public repositories and data already available (legacy data) is likely only going to be made compliant with PaCTS through dedicated synthesis efforts. Below we list some of the main issues that we encountered 465 during data synthesis. Our aim with mentioning these is to raise awareness of how the lack of standardisation affects data synthesis and thereby to encourage best practice in data reporting. We encourage researchers and also reviewers to treat data handling not as an afterthought, but as an integral part of their study. After all, compared to generating the data, data handling is not a timeconsuming task. Time spent on proper documenting and archiving is not wasted as it facilitates reuse of 470 data and enables scientific progress in our field.
Disambiguate core names: An apparently trivial, but surprisingly common, first-order issue is that core names are inconsistently archived. Different names for the same core arise from differences in hyphenation, truncation of (long) names, minor variations of the same name that can, with expert 475 knowledge, be linked, but also the use of altogether different names for the same core, e.g. reflecting differences in the labelling during an expedition and in the repository. This naming confusion renders it difficult to combine data sets from the same core, especially in an automated way, and to assess the uniqueness of time series from the same core for dereplication. We recommend to use the full name as indicated in the cruise report where the core was first described. 480 Standardise vocabularies: Even though a vast amount of palaeoclimate data is available in public repositories, the lack of standardisation of parameter vocabularies hinders efficient data processing. This problem is clearly illustrated by the fact that for this synthesis long synonym lists needed to be generated in order to group parameters. Each parameter in this synthesis had 10s-100s of different names, of which 485 many were unique (Fig. 7). This issue can be partly addressed by a consistent separation of parameter and attribute names (e.g. parameter: δ 18 O, species: G. bulloides, instead of a single parameter "d18OGbul"), but even that calls for a standardisation of parameters and attributes. Report sampling depth: All data reported here are based on measurements on discrete samples from a 490 specific depth interval in a given core. Therefore, the synthesis requires information on the position of each sample. Since ages of the samples are always estimates and may differ among studies of the same archive, unambiguous information on sample depth is essential to reproduce and update the time series.
Despite this, many studies fail to report sample depth and instead report only age. This problem is worse for spliced records that rely on a composite depth scale. Splicing approaches are often opaque and 495 original sample depths, or sample codes that identify unique samples, are not always available. We recommend therefore that sample depth should be essential for palaeodata time series from (marine) sediments and that sample labels are archived for spliced records. This includes (I)OPD/DSDP sample labels (in full) or IGSN (if available).

500
Publish raw data: To ensure reproducibility of inferred parameters (in this synthesis only temperature) and to ensure harmonisation or updating of the calibration, the raw measured data are needed. Provided that the calibration is known, measured data can in some cases be calculated from the inferred data, but this is impossible for data based on microfossil transfer functions. Related to this is unclear information about the calibration that was used, particularly if the publication describing the calibration includes 505 multiple different equations.

Include metadata:
To assess ecological imprints on proxy signals (recording bias), temperature estimates as well as oxygen and carbon isotope data based on planktonic foraminifera also require that key metadata, such as species name are archived in a standardised way. This is not universally done and for 510 example the species information is often only available in the original publication. Additional information to assess the uncertainty of proxy measurements, such as foraminifera shell size, the sample size (e.g. number of shells, concentration of alkenones), reproducibility of repeat measurements was often not available from the paper nor from the archived data and we encourage the archiving of such data in a standardised way. A similar issue applies to chronological data. Thanks to a longer history of reporting 515 standards, radiocarbon (Stuiver and Polach, 1977) (meta)data is often rather complete. However, this information is often not included alongside the digitally available data and for this synthesis a large part of the radiocarbon (meta)data had to be scraped from literature. In our age modelling approach, we took reservoir age uncertainty into account, using data derived from the modelled reservoir ages (see section 2.3). Alternative approaches are hindered by the fact that reservoir age uncertainty is almost never 520 reported.
Avoid redundancy: A considerable amount of time was spent on dereplication of time series of the same parameter from the same core that were archived multiple times. Repeat archiving happens when data are reused or, less commonly, updated. The dereplication task is not made easier by (incomplete or 525 inconsistent) metadata reporting and can be avoided through better linking of existing data sets, when instead of re-uploading the data the DOI/URL of the original data is provided.
Help rescue dark data: This synthesis is based on data that are publicly available, yet many/some palaeoceanographic time series are not archived in public repositories. Even though the proportion of 530 this so-called dark data is unknown, they likely affect every branch of palaeoceanography and as a result palaeoceanographical data syntheses cannot be exhaustive. This problem is clearly exacerbated for syntheses relying on automated data mining. There are several shades of dark data, each requiring their own approach to retrieve and make them available. Some data is only partially available, for instance data sets that lack sample depths. Such data only require additional data to make them re-useable. This 535 additional data can sometimes be calculated, or obtained from cross-referencing different data files, but in many cases will need to be retrieved from the data producers. Other data sets, in particular those from before the digital age, are presented in tables in the original publications. Progress has been made with digitising those data sets (especially PANGAEA), but this work is not finished and more effort is needed to make this data available to the community. There is also data that is used in publications, but is not made 540 available in any way (print or digital). This third shade of data can so far only be obtained from the original data producers or authors of the original publication, or if this proves impossible, needs to be digitised from graphs. Digitisation inevitably leads to loss of accuracy of the data and a data set retrieved in this way should be flagged. A final category of dark data consists of data that is not part of a publication. Such data sets can be made publicly available and be associated with a DOI to ensure 545 traceability. To reward data sharing the use of data citations needs to be encouraged and data citations should be included in the evaluation of a researcher's impact.

Acknowledgments
This research was supported by PALMOD, the German palaeoclimate modelling initiative (www.palmod.de). PALMOD is part of the Research for Sustainability initiative (FONA; www.fona.de) 550 funded by the German Federal Ministry of Education and Research (BMBF). We thank colleagues in working group 3 from PALMOD for discussion on data and metadata requirements and Rangelys Sorrentino for help with data processing. We also thank Blanca Ausin and an anonymous reviewer for their constructive feedback on an earlier version of this manuscript.

Data availability statement
An overview of all datasets used in this synthesis, including URLs to the data can be found in Table 5 (https://seafile.zfn.uni-bremen.de/f/a5732554bc91413780a3/). The PALMOD 130k marine palaeoclimate data product can be downloaded in R, LiPD and netCDF format at https://doi.pangaea.de/10.1594/PANGAEA.908831. The data can also be visualised and downloaded in 560 LiPD and csv formats at http://lipdverse.org/PalMod/current_version/.      Figure 1: Data synthesis approaches. In the expansion approach the database size increases slowly as records are added. Database size follows an opposite pathway using reduction approach and reaches a stable size quicker, with less effort. Since the expansion approach is not restricted to data that is available in the public domain, this approach may lead to a database that includes data that is not publicly available (dark data). The reduction approach on the other hand is arguably more objective, can be 590 automated and is therefore more efficient. This approach also encourages good data stewardship.  Cibicidoides. For all parameters the median resolution is more than 1 data point per 1,000 years.  Figure 6: Average temporal spacing between chronology tie points (AvgPoints). Tie points are also split in 14 C, absolute (tephra/palaeomagnetics) and tuned. The majority of the time series in the synthesis have a spacing of chronology tie points that is below 5,000 years.