Late Quaternary Climate Variability and Change from Aotearoa New Zealand Speleothems: Progress in Age Modelling, Oxygen Isotope Master Record Construction and Proxy-Model Comparisons

: We re-evaluated speleothem isotope series from Aotearoa New Zealand that were recently contributed to the Speleothem Isotopes Synthesis and AnaLysis (SISAL) database. COnstructing Proxy Records from Age Models (COPRA) software was used to produce Bayesian age models for those speleothems. The new age modelling helped us examine Late Quaternary temporal coverage for the national speleothem network, and also supported our exploration of three di ﬀ erent isotope master record generation techniques using Holocene δ 18 O data from Waitomo. We then applied the output from one of the isotope master record techniques to test an application case of how climate transfer functions can be developed using climate model simulated temperatures. Our results suggest Holocene δ 18 O trends at Waitomo capture air temperature variations weighted toward the primary season of soil moisture (and epikarst) recharge during winter. This interpretation is consistent with the latest monitoring data from the Waitomo region. Holocene δ 18 O millennial-scale trends and centennial-scale variability at Waitomo likely reﬂect atmospheric circulation patterns that concomitantly vary with surface water temperature and the isotopic composition of the Tasman Sea. A climate model simulation context for the Holocene millennial-scale trends in the Waitomo δ 18 O isotope master record suggest that site is sensitive to changes in the subtropical front (STF) and the Tasman Front. Our comparison of isotope master record techniques using Waitomo δ 18 O data indicate that caution is needed prior to merging δ 18 O data series from di ﬀ erent caves in order to avoid time series artefacts. Future work should incorporate more high-resolution cave monitoring and climate calibration studies, and develop new speleothem data from northern and eastern regions of the country. Previous work on New Zealand speleothems has deﬁned environmental shifts and climate transitions for the last glacial-interglacial cycle. However, a national-scale speleothem evaluation has not yet been undertaken using Bayesian age modelling principles. This study re-evaluates New Zealand speleothem stable isotope series that have recently been generated for SISAL (Figure and Table 1) using COnstructing Proxy Records from Age Models (COPRA) software [31], which has similar capabilities to other Bayesian approaches (e.g., BACON, StalAge). Replication of speleothem isotope signatures within individual caves, between caves located in close geographic proximity, and between karst blocks in di ﬀ erent New Zealand regions remains challenging. Previous e ﬀ orts have focused on combining isotopic data into “isotope master records” in an attempt to follow commonplace practices that replicate patterns at multiple spatial levels, which can increase conﬁdence in interpreting environmental change signals [32]. In principle, replication provided from speleothem isotope master records can diminish idiosyncrasies and edaphic factor inﬂuences on cave system drip water supplying any one speleothem. This holds true even for samples that are located within close proximity [33], as demonstrated recently with dripwater isotopes at Waipuna Cave, Waitomo [34]. Therefore, one potential beneﬁt of an isotope master record is to enhance common signals between records that are of climatic origin.


Abstract:
We re-evaluated speleothem isotope series from Aotearoa New Zealand that were recently contributed to the Speleothem Isotopes Synthesis and AnaLysis (SISAL) database. COnstructing Proxy Records from Age Models (COPRA) software was used to produce Bayesian age models for those speleothems. The new age modelling helped us examine Late Quaternary temporal coverage for the national speleothem network, and also supported our exploration of three different isotope master record generation techniques using Holocene δ 18 O data from Waitomo. We then applied the output from one of the isotope master record techniques to test an application case of how climate transfer functions can be developed using climate model simulated temperatures. Our results suggest Holocene δ 18 O trends at Waitomo capture air temperature variations weighted toward the primary season of soil moisture (and epikarst) recharge during winter. This interpretation is consistent with the latest monitoring data from the Waitomo region. Holocene δ 18 O millennial-scale trends and centennial-scale variability at Waitomo likely reflect atmospheric circulation patterns that concomitantly vary with surface water temperature and the isotopic composition of the Tasman Sea. A climate model simulation context for the Holocene millennial-scale trends in the Waitomo δ 18 O isotope master record suggest that site is sensitive to changes in the subtropical front (STF) and the Tasman Front. Our comparison of isotope master record techniques using Waitomo δ 18 O data indicate that caution is needed prior to merging δ 18 O data series from different caves in order to avoid time series artefacts. Future work should incorporate more high-resolution cave monitoring and climate calibration studies, and develop new speleothem data from northern and eastern regions of the country.

Focus of This Study
New Zealand provides critical Southern Hemisphere palaeoclimate data and perspectives about Late Quaternary environmental change [1,2]. Atmospheric circulation regimes that influence this region produce strong, heterogeneous regional hydroclimate impacts from key modes of variability like the El Niño-Southern Oscillation, Madden-Julian Oscillation, the Interdecadal Pacific Oscillation, and the Southern Annular Mode [3][4][5][6][7]. As such, understanding pre-instrumental climate history from New Zealand speleothems can help contextualize global change patterns and clarify teleconnections that operate across a range of spatiotemporal scales [8].
In this study, we provide an overview of New Zealand speleothem data that our team has provided to the Past Global Changes (PAGES) Speleothem Isotopes Synthesis and AnaLysis (SISAL) initiative [9] ( Table 1). New Zealand speleothem data (see Table 1) comprise a significant proportion of the Southern Hemisphere mid-latitude data contribution, and fewer than six caves outside of New Zealand located >30 • S were initially included in the SISALv1 database synthesis [9]. Despite speleothem investigations beginning in the 1960s [10] with extended records produced since the late 1990s [11,12], there are only a handful of published speleothems from each of the main New Zealand karst areas. New Zealand speleothems are expected to have different climatic signals from tropical Australia [13][14][15] and south Australia [16] speleothems, so their contribution toward understanding regional climate variability and long-term environmental conditions are important. Moreover, New Zealand hosts a wide range of natural archives with proxy-based reconstructions that can be compared to speleothem isotope records, and therefore they can provide a useful counterpart to Northern Hemisphere perspectives of Quaternary change [17][18][19]. Table 1. A summary of records for New Zealand cave sites reviewed for this study. Values for SISAL site id and entity ID entries for data series that are not included in either the SISALv1 [9] or SISALv2 [20] database are empty. Min year and Max year are mean age model results from COPRA that show the youngest and oldest age (rounded) of isotopic data pairs (δ 13 C & δ 18 O) for each time series. All New Zealand data series that are presented here passed screening protocols used for data accepted in the SISALv1 and v2 databases. BP = before 1950; all speleothem age models were aligned relative to this temporal datum despite being analysed in different years. See Figure 1 for the locations of SISALv1 and v2 speleothem data contributions and Figure 2 for definition of New Zealand climate regions. NA = Not yet available. NS = not suitable for palaeoclimate study.   [35]). Purple circles indicate sites that were included in SISALv1 [9], and green triangles indicate sites that were included in the SISALv2 database [20]. More details for New Zealand regional setting and climate conditions are included in Figure 2.

Site
Using the Waitomo δ 18 O records recently gathered for SISAL, we demonstrate the utility of COPRA Monte Carlo age model simulations to explore some isotope master record generation techniques. These techniques may be applicable for sites where multiple speleothem isotope series  [35]). Purple circles indicate sites that were included in SISALv1 [9], and green triangles indicate sites that were included in the SISALv2 database [20]. More details for New Zealand regional setting and climate conditions are included in Figure 2.
Using the Waitomo δ 18 O records recently gathered for SISAL, we demonstrate the utility of COPRA Monte Carlo age model simulations to explore some isotope master record generation techniques. These techniques may be applicable for sites where multiple speleothem isotope series are located in close proximity. A new Waitomo master δ 18 O isotope master record spanning the Holocene that we generated was then evaluated in a direct comparison with climate model simulations in order to surmount a situation where modern climate calibration data are limited. Our discussion focuses on the pros and cons of the isotope master record techniques we have tested and contextualises changes in the Waitomo δ 18 O isotope master record using terrestrial, marine and climate model simulation data. We conclude by suggesting future speleothem research directions, and indicate spatial gaps to target to improve national-scale speleothem coverage within and beyond the Holocene.

Physical Geography of New Zealand
Aotearoa New Zealand (referred to here after as New Zealand) is an island archipelago located in the southern mid latitudes between 35 • S and 47 • S to the southeast of Australia and adjacent to the Tasman Sea (Figures 1 and 2). Permanent human presence in this region is limited to the last 800 years [36], so many long-term New Zealand palaeoarchives were formed in the absence of anthropogenic landscape modification. Karst terrain is distributed across the length of both main islands, and confined in narrow bands that are orientated parallel to the strike of prominent northeast-southwest axial ranges [37,38].
The main axial mountain ranges intersect the prevailing westerly mid-latitude flow [39], which lends to strong west-east precipitation gradients and pronounced orographic influences on regional climates [40] that encapsulate a wide range of environments (including subtropical to glacial). The latitudinal range of both main islands means tropical to polar weather influences impact the country throughout the year [6,[41][42][43][44][45]. Teleconnections to several modes of climate variability results in highly changeable synoptic systems [40,45] that inter-and intra-annually alter temperatures, rainfall patterns and extreme environmental states (e.g., soil moisture deficit, drought, etc.).
Regional ocean circulation is complex around New Zealand [18] (Figure 1). Warm subtropical waters flow into the region at~32 • S along the Tasman Front, which is sourced from the warm waters of the East Australian Current (EAC). The EAC feeds the East Auckland Current (EAuC), which flows around the northern edge of the North Island. As the EAuC flows around the eastern tip of the east coast it forms the East Cape Current, which flows south along the east coast. The Tasman Sea to the west of the North Island is made up of Tasman Sea Central Water (TSCW). The TSCW is a combination of subtropical waters that flow to the east in the Tasman Front, to the south as part of the East Australian Current extension, and the waters associated with the broad subtropical front (STF) in the south Tasman Sea. The South Island of New Zealand is influenced by the STF, and also by the cool subantarctic waters of the Southland Current that flows north along the east coast shelf. This means New Zealand's land masses are influenced by both Subtropical and Subantarctic waters. lends to strong west-east precipitation gradients and pronounced orographic influences on regional climates [40] that encapsulate a wide range of environments (including subtropical to glacial). The latitudinal range of both main islands means tropical to polar weather influences impact the country throughout the year [6,[41][42][43][44][45]. Teleconnections to several modes of climate variability results in highly changeable synoptic systems [40,45] that inter-and intra-annually alter temperatures, rainfall patterns and extreme environmental states (e.g., soil moisture deficit, drought, etc.). , while main view shows mean annual precipitation and the distribution of karst terrain, locations where speleothems have been gathered, and boundaries for homogenous regional climate districts (following the analysis of [40,41]); (right) mean annual temperature and locations of marine records located proximal to land. 1981-2010 climatic averages are shown for terrestrial climate courtesy of National Institute of Water and Atmospheric Research Ltd. Climate districts: NNI-Northern North Island, ENI-Eastern North Island, SWNI-Southwest North Island, NSI-Northern South Island, WSI-Western South Island, ESI-Eastern South Island. Grey dashed lines represent simplified spatial patterns for δ 18 O of ocean surface water based on gridded interpolated data [46]. place in Australasia east of Australia (AUS), while main view shows mean annual precipitation and the distribution of karst terrain, locations where speleothems have been gathered, and boundaries for homogenous regional climate districts (following the analysis of [40,41]); (right) mean annual temperature and locations of marine records located proximal to land. 1981-2010 climatic averages are shown for terrestrial climate courtesy of National Institute of Water and Atmospheric Research Ltd. Climate districts: NNI-Northern North Island, ENI-Eastern North Island, SWNI-Southwest North Island, NSI-Northern South Island, WSI-Western South Island, ESI-Eastern South Island. Grey dashed lines represent simplified spatial patterns for δ 18 O of ocean surface water based on gridded interpolated data [46].
In summary, the combination of geographic position and physical geography makes New Zealand a useful location for investigating atmospheric and oceanic circulation changes using Quaternary speleothem archives [22,24,47,48].

New Zealand Records in SISAL and Additional Data
The potential of New Zealand speleothem records to extend our understanding of past climate has been demonstrated for more than a half a century [10]. Since that time, several New Zealand cave records have provided the basis for glaciation timing [11], evaluating regional environmental variability [12,24], framing climate event stratigraphy integrations (e.g., [49]), testing inter-hemispheric change mechanisms [17,22], and reconstructing atmospheric circulation patterns [21,25,27,40,50,51]. These previous studies have also helped to improve the theoretical framework for interpreting New Zealand speleothem stable isotope records ( Figure 3). has been demonstrated for more than a half a century [10]. Since that time, several New Zealand cave records have provided the basis for glaciation timing [11], evaluating regional environmental variability [12,24], framing climate event stratigraphy integrations (e.g., [49]), testing interhemispheric change mechanisms [17,22], and reconstructing atmospheric circulation patterns [21,25,27,40,50,51]. These previous studies have also helped to improve the theoretical framework for interpreting New Zealand speleothem stable isotope records ( Figure 3).  [24]). The influence of the C3/C4 ratio over interglacial/glacial timescales * is expected to be nominal for New Zealand speleothems because there is only one reported C4 plant that is a coastal sand dune specialist and is not found over cave sites. See Appendix B for details about climatic conditions and environmental effects that influence δ 18 O and δ 13 C in New Zealand speleothems.
For SISAL, New Zealand speleothem data generated from the late 1960s onward were reviewed (see references in Table 1) and submitted to the working database in multiple steps. A total of 40  [24]). The influence of the C 3 /C 4 ratio over interglacial/glacial timescales * is expected to be nominal for New Zealand speleothems because there is only one reported C 4 plant that is a coastal sand dune specialist and is not found over cave sites. See Appendix B for details about climatic conditions and environmental effects that influence δ 18 O and δ 13 C in New Zealand speleothems.
For SISAL, New Zealand speleothem data generated from the late 1960s onward were reviewed (see references in Table 1) and submitted to the working database in multiple steps. A total of 40 speleothems were gathered and screened, and only data series with two or more geochronology dates were included in this study ( Figure 1). A total of 29 unique speleothems from 17 caves with δ 18 O and δ 13 C records were found to be suitable for further analysis. The locations of these speleothem records are shown in Figures 1 and 2, and temporal details for each speleothem are summarised in Table 1 and Figure 4. Information about speleothems that were considered unsuitable for palaeoclimatology using stable isotopes, that do not have geochronology data, or are not yet available for palaeoclimate analysis are also listed in Table 1. Further details about the regions the speleothems were collected from are reiterated from previous work [25] in Appendix A. δ 13 C records were found to be suitable for further analysis. The locations of these speleothem records are shown in Figures 1 and 2, and temporal details for each speleothem are summarised in Table 1 and Figure 4. Information about speleothems that were considered unsuitable for palaeoclimatology using stable isotopes, that do not have geochronology data, or are not yet available for palaeoclimate analysis are also listed in Table 1. Further details about the regions the speleothems were collected from are reiterated from previous work [25] in Appendix A. . Temporal coverage of New Zealand speleothem records submitted and/or prepared for SISAL (covering MIS4-present) during this study. Note, some records are not currently included in the SISAL database but remain in preparation (see details in Table 1). Shading denotes the average temporal resolution for consecutive isotopic samples in years. Groupings of speleothem entities by cave location are ordered top to bottom by latitude. Figure 4. Temporal coverage of New Zealand speleothem records submitted and/or prepared for SISAL (covering MIS4-present) during this study. Note, some records are not currently included in the SISAL database but remain in preparation (see details in Table 1). Shading denotes the average temporal resolution for consecutive isotopic samples in years. Groupings of speleothem entities by cave location are ordered top to bottom by latitude.

Dating, Age Models and Isotopic Analyses
Use of radiocarbon ( 14 C) to produce chronologies for New Zealand speleothems was commonplace through the 1990s [12]. In addition, 210 Pb analyses were also employed on soda straw calcite formations to attempt high-resolution isotopic data calibrations for palaeoclimate reconstructions [12]. Geochronology of speleothems in this study were provided via Uranium-Thorium (U-Th) disequilibrium dating following previously outlined methods [21,23,24,52]. All speleothem isotope data that were generated at NIWA (see Table 1, [21,22,24,25], and Williams unpublished) had associated U-Th analyses generated using thermal ionisation mass spectrometry (TIMS) at the University of Queensland [53]. Geochronology analyses on speleothem data generated by University of Melbourne (ED1 and MD3) were established using methods previously described [23,30]. Mean and median errors for all data used in the age modelling is approximately 1.0% and 0.5%, respectively (N (the number of samples) = 148; std. dev. = 1.4%), and most of the U-Th dates that have errors >2% come from Waitomo and Hawke's Bay.
Both 14 C and U-Th approaches have inherent issues that are recognised for some New Zealand caves [28]. However, robust numbers of dates in each speleothem helps to minimise this issue, especially in the context of Bayesian age modelling. For this study, individual speleothem records were age modelled in COPRA software (version 1.15) inclusive of the depth data, geochronology constraints (U-Th ages with errors) and stable isotope values [31]. This effort has updated previous work that did not consistently publish full age modelling data with stable isotope series. All speleothem time series (except RK05-3 [27]) were age modelled using the Piecewise Cubic Hermite Interpolation Polynomial (PCHIP) and 1000 Monte Carlo simulations to generate mean and median chronologies with confidence intervals for each isotopic data point. The 1000 age model iterations were also retained from COPRA for use in master isotope construction approaches (see Section 2.4).
All of the oxygen isotopic data generated at NIWA and University of Melbourne (the majority of the data series) used powdered samples of finely milled speleothem CaCO 3 run on Finnigan MAT251 and MAT252 devices equipped with a Kiel automated individual-carbonate reaction chamber (which were combined with H 3 PO 4 to derive CO 2 ). All speleothems examined in this study were tested in the original studies for isotopic equilibrium via analysis of the axial isotopic profile. While we recognise that these tests are not definitive [54,55], interpretations of isotope records (in particular from Waitomo) are lent greater veracity by the general agreement of multiple speleothem records within a single regional setting.

Controls on Speleothem Growth and Interpretation of Speleothem Stable Isotope Signals
The supply of water to karst bedrock (either regularly or intermittently) is requisite for speleothem formation [55,56]. Across New Zealand, rainfall events can be frequent for most regions year-round [41] despite some prolonged, intra-seasonal dry spells. New Zealand speleothem calcite δ 18 O and δ 13 C can be influenced by a range of processes-some of which operate on microscopic scales, and others that arise from broad environmental influences dictated by global climate variability and change [12,21,24]. Key principles about the drivers of New Zealand speleothem δ 18 O and δ 13 C variability on centennial to orbital time scales, and details about interpreting δ 18 O and δ 13 C as palaeotemperature and palaeohydrology proxies, are based on previous work [12,21,22,24,25,30,51]. We provide a summary diagram below ( Figure 3) and additional information from previous work in Appendix B.
Briefly, the primary controls on precipitation δ 18 O values in New Zealand include temperature [57], moisture sources [58], and rainfall amount [34] (see more details in Appendix B). Data from a network of precipitation isotope monitoring sites in New Zealand demonstrate a strong control of temperature, altitude, and precipitation amount on precipitation δ 18 O values [57]. Highest mean annual δ 18 O values are found at lower latitudes and altitudes and near the coasts, with lowest δ 18 O values in the lee of the Southern Alps on the South Island where northwesterly air masses experience the greatest fractional rainout [59]. In the Waikato area near the Waitomo cave sites, seasonal precipitation δ 18 O values are inversely correlated with rainfall amount, which is consistent with the "amount effect" [34]. In the Waitomo area, cave infiltration δ 18 O values are similar to mean annual precipitation δ 18 O values [34]. Waitomo speleothems may therefore be plausibly interpreted to reflect changes in precipitation δ 18 O, which can vary due to changes in temperature, rainfall amount, or moisture source on decadal and longer time scales.
But what are the primary long-term controls on precipitation δ 18 O changes that are needed to best interpret Holocene speleothem δ 18 O variations? One approach to aid speleothem δ 18 O interpretation is to establish a calibration based on modern climate and precipitation δ 18 O variations. This approach is difficult to achieve in New Zealand because the epikarst overlying the cave smooths the δ 18 O signal of drip waters relative to above-cave precipitation [60] (and in some cases by an undetermined amount). It may also be unclear which statistical fit between precipitation δ 18 O and climate variables is most appropriate. Therefore, resulting calibrations are commonly limited to short time periods of precipitation δ 18 O observations, and stationarity of the climate/δ 18 O relationship often has to be assumed. Alternatively, calibration of speleothem δ 18 O values to modern climate observations can provide a transfer function (e.g., [61]). However, the speleothems in this study either grew too slowly or were sampled too coarsely to yield a valid calibration. We have therefore developed a calibration between the isotope master speleothem δ 18 O time series using modelled climate variables. This approach bypasses the uncertainties of applying calibrations from the limited instrumental observation period to the speleothem data.
The primary factors that control calcite δ 18 O (Figure 3; further details in Appendix B) indicate that uncorrected δ 18 O isotopic data for many New Zealand speleothems that extend beyond the mid-Holocene contain a global ice volume signature [22,47]. Continental ice sheet expansion and corresponding sea level lowering resulted in isotopic enrichment of δ 18 O in sea water (δ 18 O sw ) which influences δ 18 O precipitation . As such, an ice volume correction needs to be applied to New Zealand δ 18 O calcite speleothem records that are older than the mid-Holocene. This correction is applied under the assumption that precipitation δ 18 O values tracked the change to δ 18 O sw , and that it is required before palaeoenvironmental interpretations are made from δ 18 O calcite . In this synthesis, we choose an existing relative sea level model [62] that includes a sea level nadir estimate of 123m below present day during the global last glacial maximum (LGM), and follow the assumption that 1.2 +/− 0.1% represented the average δ 18 O sw change over the last four terminations [63]. Isotopic enrichment of δ 18 O calcite was also assumed to relate linearly to sea level change. An extant sea level record [62] was linearly interpolated using those assumptions to derive an annually-resolved ice-volume δ 18 O correction factor, which was then applied to uncorrected isotope master record δ 18 O to derive a corrected series, denoted δ 18 O ivc .

Isotope Master Record Approaches
We used COPRA age model data and multiple approaches to demonstrate how master isotope records can be constructed, with a specific focus on Holocene δ 18 O trends and variability at Waitomo Caves. Waitomo has many caverns in close proximity that have yielded speleothems used in palaeoclimate reconstructions [12,24]. We focused on the Waitomo data specifically due to the abundance of individual speleothems that were available for isotopic data comparisons (and other reasons discussed below). Putting microscale conditions aside, we treat the climate of the Waitomo region as a homogenous entity, because it is a small area and therefore we expect isotopic trends between speleothem records will be similar [64]. Previous work demonstrates the potential to replicate regional palaeoclimate signals [12,24] and emphasise common patterns using isotope master records developed from Waitomo speleothem isotopic series [47,48]. However, little work has been reported on comparing different methods that combine isotopic series into master records, which we address here.
First, median age models derived from COPRA were applied to individual speleothem δ 18 O time series from the COPRA output, which accounts for age model uncertainty. Then, all isotopic data points from all series were interleaved in chronological order, plotted together, and fitted with a 5-point running mean. This approach was essentially used by Williams et al. [22,24,47]. The depiction of this type of master series, however, does not distinguish the independent age uncertainty estimates related to individual isotopic data points from each speleothem during the process of combining all the data into one record. Second, we used 2.5th and 97.5th percentile uncertainty bands for the COPRA age output that was derived via a Monte Carlo simulation and 1000 age models produced for each speleothem isotopic series to create a "binned" isotope master record. Original δ 18 O values with a modelled age that fell within a 100-year time span were collated and averaged for individual speleothems using arithmetic mean and equal weighting, independent of isotopic value distance from central bin age. 100-year bin averages were iteratively calculated in 50-year overlapping steps for all of the individual speleothem δ 18 O series. Discrete 100-year binned δ 18 O values for all individual speleothems were then adjusted via normalisation to the Gardner's Gut 2 (GG2) speleothem series, and then all data were averaged into an isotope master record. As such, this type of isotope master record inherently includes consideration of the dating uncertainty for each speleothem.
Last, we used a Monte Carlo Empirical Orthogonal Function (MCEOF) approach to evaluate common low-frequency millennial-scale trends in coeval speleothem records. The MCEOF uses Monte Carlo-derived age model simulations constrained by geochronology uncertainties with Empirical Orthogonal Function analysis. It generates a mean time series that captures maximum common variance for independent speleothems with different age models that overlap in time. We also consider that the MCEOF approach could be used as a mechanism to merge speleothem isotope records, and we evaluated whether mean isotopic value shifts can occur with the introduction and omission of individual speleothem isotopic data series through time during master record construction. We harnessed the 1000 Monte Carlo age models produced in COPRA that were used to estimate the chronologic uncertainty of each isotopic series (see "binned" approach above). The results from the Monte Carlo age modelling for each speleothem and the common period of overlap between different speleothems also dictated a minimum time coverage for each MCEOF series. This approach was applied using all 1000 age model iterations between coeval speleothems, including as few as two and as many as five speleothem series at one time to create a MCEOF master record. During this process, each resulting EOF was rescaled to be compatible with original δ 18 O values common to each speleothem archive. The MCEOF technique, like the second one described above, includes full considerations of dating uncertainties for the isotopic values in each speleothem.

Time Series Coverage Limitations for Isotope Master Records
Multi-millennial overlaps for speleothem coverage guided pre-screening of samples for isotope master record generation. The region with the most samples is Waitomo, followed by North Westland; the fewest come from Paturau and Hawke's Bay. As such, we focused on Waitomo and North Westland as initial candidate regions to evaluate isotope master record techniques because they have the greatest data coverage.
As a data screening mechanism, we first examined δ 13 C and δ 18 O correlations from Waitomo and North Westland speleothems. We relied on first-order principles of interpreting those data conjointly ( Figure 3; Appendix B) to assess whether the oxygen isotopic signal may be strongly impacted by local water balance (the flow of water into and out of the epikarst/karst, and by way of association with effective precipitation). This assumption relies on accepting δ 13 C is a proxy for water balance (and effective precipitation) following previous work (see Appendix B). We observe a lack of correlation between δ 13 C and δ 18 O for Waitomo series (n = 7; mean r = −0.07; max r = −0.2; min r = 0.0; no significance above the 95% confidence level), suggesting the δ 18 O signal in speleothems gathered from that region may only be weakly influenced by water balance changes (see Table 2).
For North Westland, the situation is very different, with a majority of speleothems showing significant δ 18 O-δ 13 C correlation coefficients > 0.3 (n = 7, mean r = 0.35; max r = 0.72, min r = 0.0; five of seven p-values < 0.01)) (see Table 2). This suggests a strong potential imprint of water balance on δ 18 O in North Westland speleothems (if previous theoretical assertions about δ 13 C hold true for this location). When North Westland and Waitomo climates are compared, the former area receives about twice as much rainfall as the latter ( Figure 2). It is therefore not surprising that a potential rainfall amount effect (and water balance surplus) might have a more significant impact on speleothem δ 18 O in North Westland where precipitation is greater. We note that there could be an important component of rainfall source region seasonality for North Westland that may also contribute to speleothem δ 18 O signatures there. For those reasons, only a master δ 18 O series from Waitomo was pursued in this study, and we express caution for drawing temperature interpretations from North Westland primary δ 18 O data.  [65]. Preliminary climate response function analyses that have compared limited modern Waitomo δ 18 O data to local instrumental data supports this assertion [25]. We attempted a transfer function approach that compared speleothem δ 18 O to modelled climate simulation temperatures. This approach also presents a way to bypass some of the complexities of interpreting the competing influences on speleothem δ 18 O and allows a consideration of the functional δ 18 O ivc response to large-scale climate change in the pre-instrumental period. We also assumed that a constant cave temperature does not affect speleothem δ 18 O via variations in the calcite-water 18 O fractionation, which should be small relative to changes in the δ 18 [69]). Annually-resolved mean monthly temperature data from the grid cell closest to Waitomo (~38 • S and~175 • E) were obtained from each model; three-month averages (e.g., June-July-August, JJA; July-August-September, JAS; August-September-October, ASO; etc.) were also generated from the monthly climate model data. Modelled temperatures were then binned into 100-year averages staggered through time every 50 years (i.e., the same way the Waitomo δ 18 O isotope master record was constructed). The 100-year binned modelled temperatures were chronologically aligned to the Waitomo δ 18 O ivc isotope master record using 1950 CE as a starting point to enable correlations with each climate simulation data series.
Linear regressions were then used to derive a transfer function that converted Waitomo δ 18 O ivc (in per mille) to modelled centennial-average JJA temperatures (austral winter). Twice the standard deviation of the average residual for the linear regression was employed as the 95% confidence bound for converting δ 18 O ivc to temperature. In order to compare δ 18 O-derived Holocene temperature changes based on different transfer functions that were developed from each unique climate model simulation (which have different absolute temperatures), each temperature-transformed isotope master record was converted further to a JJA temperature anomaly relative to the last 4 ka (Late Holocene) average.

New Zealand Speleothem Coverage and Further Focus on the Holocene (MIS1;~11 ka-Present)
For New Zealand speleothem isotopic data series that were reviewed and age modelled for SISAL, only 10 records extended beyond the late glacial, and half of those limited records extended prior to the LGM (Table 1 and Figure 4). Holocene speleothem coverage includes series from both main islands (six series from Waitomo, three from Hawke's Bay, one from Paturau, three from Kahurangi, three from North Westland and three from Fiordland-Southland.) Each of these regions has a minimum of three speleothems that cover all or part of the Holocene, with sample resolution widely varying from sub-decadal to multi-centennial scales ( Figure 4). Many of these speleothems were previously used in palaeoclimate reconstructions for the Late Holocene (4 ka-present; [21]), the mid-Holocene (~7-5 ka), the last 2000 years [25], and the last glacial-interglacial cycle [38]. In this analysis using COPRA age modelling, temporal coverage for isotopic data spanning the Holocene (~11 ka-present, with all data binned at centennial resolution) ranges from a high of 100% (i.e., all speleothems covered an interval) to a low of 17% (mean = 46%; STD = 21%).
Our choice to focus in more detail on Holocene speleothem data is due to greater similarity of intra-cave δ 18 O isotopic signatures, in addition to greater data depth and better spatial density of regional coverage during the Holocene. MIS2-MIS4 data, as well as more detailed interpretations of δ 13 C, will be described elsewhere. All dates reported below are equivalent to calendar years before present (1950). We note that the isotopic profiles of speleothems collected from different locations in each cave can differ on sub-millennial time scales, possibly from different sensitivities of drip waters to surface δ 18 O variations, kinetic isotope effects, growth rate and sampling resolution differences, and/or age model uncertainties. Nevertheless, the millennial-scale δ 18 O signals observed for some closely spaced Holocene records (including the Waitomo speleothem data) appear robust ( Figure 5).

Speleothem Isotope Master Record Techniques Using Holocene δ 18 O Data from Waitomo Caves
Inter-speleothem correlations for the 100-year binned Waitomo δ 18 O data were derived for all overlapping series using a matrix containing two sets of 10 3 COPRA age modelled results. This data pre-screening approach produced 10 6 inter-series correlations that indicate most of the δ 18 O records have a common signal (see Table 2). However, we rejected Max's Cave from further consideration in constructing a master δ 18 O record due to poor overall median inter-series correlations (r = 0.19). Rejection of Max's Cave δ 18 O is one difference from previous work for Waitomo [12,24,47], which had included those data in a isotope master record.
Several methods were then tested to combine speleothems into an isotope master record, which produced subtly different results ( Figure 5). Use of the median age model from COPRA and a 100-year bin technique produced time series with variability and structure that is similar to the 5-point running mean technique of Williams et al. [24,47,48]. We note, however, that some 100-year time steps have no data coverage due to the constraints arising from dating uncertainties and sampling resolution incorporated within the COPRA age modelling.  The MCEOF technique that used temporally-overlapping COPRA simulated age models produced more emphasis on common low-frequency signals between speleothems (see Figure 5; >85% variance explained in each MCEOF iteration). MCEOF isotope master record coverage is always limited to the period of common overlap for the series incorporated into each master record. Because there is a blend of 1000 iterations of each series into the master, the result also has much smoother trends and reduced variability relative to the original data. The MCEOF technique also shows the influence of how individual isotope series affect the average δ 18 O value in an isotope master record ( Figure 5, third panel).
The unadjusted Waitomo δ 18 O isotope master records that span the mid-to-late Holocene (e.g., 8 ka-present) show millennial-scale trends that are confined well within a~1% span. The 5-point running mean and 100-year binned records retain centennial structure, typically on the order of up tõ 0.2 to 0.5% fluctuations from the local mean that show episodic excursions and, occasionally, more or less pronounced variability. More positive mean δ 18 O values are generally observed for most of the early Holocene (10.8-6.5 ka) and part of the Common Era (last 2 ka; including a Late Holocene peak between 1.2 and 0.8 ka) for all of the Waitomo isotope master records. This assertion holds even when the Waitomo δ 18 O isotope master records are adjusted for ice volume effects (see Figure 6). The most consistently negative isotopic levels across millennial scales, representative of the mean climate state, are observed at~3 ka. There is less isotopic variability across most of the early-to-mid Holocene, but this is potentially an artefact of fewer speleothem isotope series and fewer data points that contribute to each isotope master record during that time. When there are an elevated number of isotopic data points that contribute to each isotope master record during the early Holocene, the variability appears equivalent to the mid-to-late Holocene ( Figure 5). In the next section, we adopt the 100-year binned isotope master record for further analysis after adjusting it for ice volume changes (and we explain further in the discussion why this record may be best to use).

Comparison of Waitomo δ 18 O Master Record to Climate Model Simulation Data
Correlations between transient climate model simulated monthly temperatures and Waitomo δ 18 O ivc master series binned at 100-year increments indicate a significant positive δ 18 O relationship during most of the austral winter months (JJA) for three out of four models, which aligns to the season of main hydrologic recharge (late autumn-winter-early spring [34,64]). JJA average modelled temperatures (austral winter) improved the strength of correlations between the climate model output and the Waitomo δ 18 O ivc master series (see Table 3). The most negative Holocene δ 18 O ivc values largely occur during the mid-to-late Holocene between 3 and 5 ka, when modelled JJA temperatures are the coolest. The least negative values occur during the early Holocene, and within the last 2 ka when modelled temperatures were warmest ( Figure 6). One of the models (MPI-ESM1.2) used to transform the Waitomo δ 18 O ivc isotope master record to palaeotemperature shows little temperature change through the Holocene. Each of the remaining climate model simulations show a different scale of low frequency change, and a distinct relationship between δ 18 O and temperature, that ranges from 0.23 • C to 0.7 • C per mille.
The common millennial-scale trend for the Waitomo δ 18 O ivc isotope master record and modelled JJA temperature ( Figure 6) shows some close associations between Holocene temperature changes at Waitomo and solar irradiance flux [70] for a portion of the Holocene (addressed in the discussion in more detail). While JJA millennial-scale temperature trends across the Holocene are prominent, there are also clear centennial-scale temperature departures (collectively indicated by the transformation of δ 18 O to temperature) of up to +/−0.3 • C relative to the millennial-scale trends. Although we interpret these centennial-scale variations of δ 18 O ivc as related to temperature, the more abrupt shifts may also be plausibly related to temperature-independent shifts linked to atmospheric circulation or changes in rainfall amount. We note that inter-annual to multi-decadal variability for JJA seasonal temperatures of up to +/−2 to 3 • C, relative to the evolving millennial-scale average, are also exhibited by the transient climate model simulations. These patterns are not expected to be observed in the Waitomo δ 18 O isotope master record because of how it was constructed, in addition to the (coarse) sampling resolution relative to the deposition rate for the speleothem materials analysed.
Time slice simulation output for 6 ka from PMIP2 [71] are more numerous than available transient Holocene simulations. They collectively show near-average or slightly below average temperatures for June-September ( Figure 6) in the North Island and the Tasman Sea west of Waitomo. The PMIP2 multi-model ensemble average and spread for 6 ka temperature anomalies (Figure 6 inset) overlaps the transient simulation results. The temperatures derived for the Waitomo δ 18 O ivc master series appear displaced slightly higher than the transient model output for 6 ka. However, the spread for the (related) PMIP2 models and the scale of associated uncertainty from the transfer function application suggest the temperature-transformed isotope master record and a range of independent models agree during this time.
be plausibly related to temperature-independent shifts linked to atmospheric circulation or changes in rainfall amount. We note that inter-annual to multi-decadal variability for JJA seasonal temperatures of up to +/−2 to 3 °C, relative to the evolving millennial-scale average, are also exhibited by the transient climate model simulations. These patterns are not expected to be observed in the Waitomo δ 18 O isotope master record because of how it was constructed, in addition to the (coarse) sampling resolution relative to the deposition rate for the speleothem materials analysed.

What Can New Zealand Speleothems Gathered for SISAL Tell Us about the Late Quaternary?
The isotopic signatures evaluated in previous New Zealand speleothem studies [47,48] have been compared with Late Quaternary vegetation and marine climate records [72][73][74][75][76][77][78][79][80]. Contemporary ecology studies that support Late Quaternary pollen interpretations establish an intimate link between New Zealand plant ecosystems and multiple climate variables [81]. During the LGM between~32 and 18 ka, previous work shows more negative speleothem δ 18 O isotopic values [24,25,49] occur alongside some of the most positive δ 13 C values. This suggests those dual proxies reflect colder and drier conditions, respectively, using theoretical understanding of speleothem cave calcite deposition ( Figure 3). LGM isotopic signatures for some New Zealand speleothems show anomalously negative δ 18 O and positive δ 13 C values. These signatures align with widespread grass/herb expansion and tall tree/thermophilous plant reduction [78,79] at a time when offshore sea surface temperatures were as much as 5 • C cooler than the Holocene [82]. There is also evidence of near-maximally expanded glaciers on the North and South Islands [83,84] during the LGM, which required a combination of both colder and drier conditions to achieve geometries that best match moraine limits [85]. Conversely, speleothem δ 18 O isotopic values equivalent to or less negative than modern levels (e.g., during the early Holocene interglacial) occur when conditions as warm or warmer than present are expected.
The correspondence between speleothem isotopes and palaeoproxy evidence over glacial-interglacial scales supports assertions that cave calcite δ 18 O and δ 13 C both have utility as environmental variability and climate change records [12,21,22,24,27,48]. In addition, New Zealand speleothem isotopic data can define climatic intervals and key transitions during the last glacial-interglacial cycle that other proxies cannot. One obvious benefit using speleothems comes from the ability to provide secure geochronology with U-Th beyond the limit of radiocarbon, when long lake/bog pollen records are difficult to date accurately. In a similar sense, long marine sediment records beyond the range of radiocarbon are also dependent on wiggle-matching oxygen isotope data from foraminifera to global oxygen isotope curves or Antarctic temperature records.

Replicated Isotopic Signals Support Creation of Isotope Master Records
Merging different isotopic data series together into master records can potentially provide a longer view than what can be achieved from one speleothem archive. In New Zealand, where significant seismic, volcanic, and glacial activity has occurred in proximity to some caves (potentially causing abrupt speleothem growth rate changes), there may be benefits with isotope master records that can surmount deposition inconsistencies and hiatuses.
The isotope master record approaches we tried are not exhaustive, but expand on previous work that applied a 5-point running mean through all available data [24,47,48]. The 100-year binned isotope master record ensures even sampling through each individual series is adhered to. Subsequently, direct comparisons between speleothems from within a karst block using this method inherently retains sample depth metadata that can be directly linked to Bayesian age modelling. The 100-year binned approach does not introduce any chronology or interpolation artefacts, and it diminishes the possibility of spurious inter-data point trend interpretations. In addition, some bins are dominated by a lone isotopic measurement. Caution is warranted for interpreting any extreme conditions in those instances; however, the retention of sample depth information using the 100-year binned approach means these situations can be clearly identified.
The MCEOF approach shows how mean isotopic value shifts can occur when certain speleothems are introduced into an isotope master record composed of multiple data series. While some previous techniques that have merged isotopic series together (e.g., 5-point running mean) probably contain minor statistical artefacts, the MCEOF approach highlights the potential for evaluating the severity of that type data inhomogeneity. Therefore, the MCEOF method has useful application for re-evaluating the validity of transfer functions and the interpretation of quantitative climate reconstructions from older master records where previously low sample depths have been recently improved. The MCEOF method also has potential application in situations when a "master" series to adjust all others to cannot be decided on. A conservatively re-scaled MCEOF isotope master record could be used to adjust overlapping isotopic data via matching isotopic means and slopes. Variance differences that arise between each adjusted series in that situation, however, would still need to be attended to using other techniques before the data are combined into an isotope master record.

Waitomo Holocene δ 18 O Isotope Master Record Temperature Context
We calibrated and transformed the 100-year binned Waitomo δ 18 O ivc isotope master record using transient climate model simulated temperatures, after recognising there were similar low-frequency millennial-scale trends in both datasets. This is not the first time that a speleothem δ 18 O series has been calibrated using palaeoclimate model information [86]. Waitomo δ 18 O ivc shows the least negative values during the early Holocene (10.8-6.5 ka) and between 1.2 and 0.8 ka. The calibration we have used, in addition to theoretical understanding, suggests these periods may have been warmer than late Holocene pre-industrial times (150-250 BP).
Early Holocene pollen records suggest tall trees were more expansive along with increased presence of frost-tender and drought-intolerant taxa, relative to the hardier species that were formerly more dominant during the LGM [81,87]. This implies the early Holocene was warm, with reduced frost and drought incidence [88] and increased humidity [73]. In addition, previous pollen work has asserted the early Holocene conditions reflect influences on the mean climate state from diminished seasonality [89,90]. The most negative Waitomo δ 18 O ivc values occur in the mid-to-late Holocene (after~4 ka) and suggest this was the interval with the coolest winter temperatures in the present interglacial. Mean δ 18 O values for the last 1000 years suggest winter temperatures at Waitomo were warm at~0.9 ka, followed by near normal conditions through the late pre-industrial period (~0.5-0.1 ka) inclusive of a small downward trend to more negative values between~0.5 and 0.2 ka. While inter-annual to multi-decadal temperature variations are not captured by our isotope master record, high resolution tree-ring data capture that type of activity [91][92][93][94][95][96], and suggest sub-centennial variability was an important climate component across the last millennium. Winter palaeotemperatures inferred from Waitomo δ 18 O isotope master record values for the last few hundred years do not appear as cold as the~2.2-5.8 ka time span (Figure 6).
New Zealand's glaciers are acutely sensitive to summer temperatures [97,98]. Negative temperature anomalies inferred from moraine-constrained glacier modelling [99,100] and alpine glacier palaeoequilibrium line altitudes [50,84,101] concur with tree-ring evidence [91] to suggest mean summer climate over the last few hundred years was cooler than present. Moraine evidence also suggests summer temperatures during significant portions of the Holocene may have been colder than the present day and the late pre-industrial average [21,84,99,[102][103][104][105][106]. The juxtaposition of two different proxy types (speleothems and glaciers) that reflect opposite seasonal sensitivities (winter vs. summer) independently reinforces assertions about reduced seasonality influences during the early Holocene (with warmer winters and cooler summers) that have been conceptually made using palynology [73,107].
Assertions about early Holocene warmth from the Waitomo δ 18 O ivc isotope master record are also supported by Holocene marine sediment proxy records of oxygen isotopes and sea surface temperatures around northern New Zealand and the Tasman Sea. Oxygen isotopes from the east coast of Australia indicate increased southward flow of the EAC from around 11 ka [108]. EAC influences on Tasman Sea temperatures during the Holocene were also likely to have been transferred by the Tasman Front to the east of New Zealand, and further propagated around the North Island via the EAC and then the East Cape Current (see Figure 1 for ocean current locations). Northeast of the North Island, SST are slightly warmer during the early Holocene and display a subtle temperature decline through the Holocene (H214 [109,110]; JP87 [111]). The δ 18 O of planktic foraminifera Globigerina bulloides for this area show lowest values between 10 and 8 ka and then higher values throughout the rest of the Holocene (H214 [109]; JP125 [110,112]). The consistent offset between the δ 18 O of G. bulloides (a mixed layer (100 m) foraminifera) and Globorotalia inflata (a deeper dwelling planktic foraminifera) indicates a well-stratified water column during the early Holocene.
Warmer temperatures also existed at~9 ka to the southeast of the North Island, followed by a steady decline of~1 • C by 4.5 ka (MD97-2121 [110,113]). The δ 18 O of G. bulloides are variable, but show a decrease around 8 ka [114], while the δ 18 O of Globigerinoides ruber (surface dwelling foraminifera) and G. inflata (which lives in the thermocline, deeper in the water column) indicate waters in this region were well-stratified throughout the Holocene [114]. Nearby, radiocarbon dates that constrain mangrove (Avicennia marina [115]) presence south of East Cape in Gisborne (south of modern limit) have also suggested warmer winter temperatures (and reduced frost) existed in the nearshore environment at the onset of the early Holocene.
West of the South Island, peak early Holocene warmth is also observed, followed by a relatively steady temperature decline, evident from a decrease in alkenones SST proxy and an increase in δ 18 O G. bulloides during the last 10 ka (SO136-11 [116]). Foraminiferal assemblages and alkenones from the South Tasman Rise (south of Tasmania) also indicate slightly warmer temperatures during the early Holocene, likely linked to the extension of the EAC and the southward shift of the STF. Cooling occurred after 8 ka, reflected in G. bulloides δ 18 O values that were achieved between 6 and 4 ka (RS147-GC7 [117]). The STF (Figure 1) was also shifted south during the early Holocene by~1 • latitude bathing the South Island of New Zealand in warmer waters [117,118], and then it rebounded to its modern position between 8 and 6 ka.

Reconciling Waitomo δ 18 O Master Record Temperature Interpretations with Marine Records, Terrestrial Records and Model Results
Continent-ocean temperature contrasts and differences between low-to-mid latitude and high-latitude Pacific and Southern Ocean temperature anomalies are obvious in both time slice and transient Holocene simulation experiments (Figure 7). They indicate New Zealand's maritime climate setting may exhibit unique signatures of change relative to other southern mid-latitude landmasses and oceanic sites. In addition, pan-oceanic spatial differences for proxy and modelled climate trends may arise, in-part, because ocean currents and marine climate inertia blend inter-seasonal temperature anomalies, and because models are coarsely resolved. Some location-specific issues, including the intricacies of ocean currents and gyres and some potential seasonal biases for proxies, can further confound direct and consistent onshore-offshore evidence comparisons for proxy-model syntheses. As such, regional-scale proxy-model comparisons are recommended [68]. The initial focus on SSTs proximal to New Zealand and regional data interpretations are therefore justified as a conservative starting point to contextualize the Waitomo δ 18 O ivc isotope master record. Below, we discuss the evidence for millennial-scale climate changes in the New Zealand region, and interpret trends for the Waitomo δ 18 O ivc isotope master record in two parts.

Did Ocean Source Region Variability Consistently Influence Holocene Waitomo δ 18 O Trends?
Early Holocene marine climate proxy evidence indicates warmer-than-modern temperatures existed around the Tasman Sea, which then cooled into the mid-Holocene (see Section 4.3). The air temperature change inferred from the model-calibrated Waitomo δ 18 O ivc isotope master record shows millennial-scale trends and centennial variability was constrained within +/−1 • C of modern conditions.
Prevailing westerly atmospheric flow that delivers rainfall to the North Island would impart Tasman Sea source water physical traits onto Waitomo's cave drip water pool held in the epikarst [34]. The early-to-mid Holocene shift toward more negative δ 18 O values at Waitomo may reflect Tasman Sea δ 18 O changes, particularly if the ocean source effect dominates over the cave temperature effect. For Waitomo speleothems, fed by drip water in a maritime climate setting where relatively unimpeded atmospheric transport pathways to the adjacent ocean upwind exist, this possibility is realistic [119]. Alternatively, the δ 18 O decrease from the early-to-mid Holocene may represent an influence of cooler temperatures on precipitation δ 18 O values, a greater contribution of high-latitude moisture masses from the Tasman Sea, or a larger precipitation amount effect. In addition, a larger component of Pacific moisture that experienced more 18 O depletion upon traverse of the North Island could explain a δ 18 O decrease; however, this mechanism is unlikely because Waitomo is located near the Tasman Sea coast.
A shift in Tasman Sea δ 18 O between the early-and mid-Holocene could have also arisen from contraction of the STF southern boundary to its present position by the mid-Holocene. However, modern ocean chemistry isotopic maps indicate a~1 • of latitude northward shift of the STF equates to only about −0.1% for δ 18 O of Tasman Sea ocean source waters (Figure 2). This suggests the early-to-mid Holocene STF rebound may have not been sufficient as the sole mechanism for influencing Waitomo δ 18 O calcite via the drip water source pool. realistic [119]. Alternatively, the δ 18 O decrease from the early-to-mid Holocene may represent an influence of cooler temperatures on precipitation δ 18 O values, a greater contribution of high-latitude moisture masses from the Tasman Sea, or a larger precipitation amount effect. In addition, a larger component of Pacific moisture that experienced more 18 O depletion upon traverse of the North Island could explain a δ 18 O decrease; however, this mechanism is unlikely because Waitomo is located near the Tasman Sea coast.  [120]). These results are from model counterparts of CCSM3-TraCE 21k, Echo-G, and ECBilt-CLIO models, respectively, noted in Figure 6.
A shift in Tasman Sea δ 18 O between the early-and mid-Holocene could have also arisen from contraction of the STF southern boundary to its present position by the mid-Holocene. However, modern ocean chemistry isotopic maps indicate a ~1° of latitude northward shift of the STF equates to only about −0.1‰ for δ 18 O of Tasman Sea ocean source waters (Figure 2). This suggests the earlyto-mid Holocene STF rebound may have not been sufficient as the sole mechanism for influencing Waitomo δ 18 Ocalcite via the drip water source pool.
Reduced EAC subtropical water inflow (warm, isotopically less negative) and weakening of Tasman Front influences between the early Holocene and mid-Holocene may be another key component to consider for Tasman Sea marine climate changes (and downwind connections to Waitomo). Modern ocean chemistry spatial relationships show a ~0.2‰ northwest-southeast δ 18 O gradient for sea water situated between Australia and the North Island. Enhanced southward EAC transport during the early Holocene, with an expanded and invigorated Tasman Front, may have pushed isotopically less negative δ 18 O waters to the south that normally are positioned to the northwest of the North Island and adjacent to Australia's east coast. In this situation, isotopically less negative ocean waters could have been advected via westerly drift into the central Tasman Sea. When the EAC and inflow along the Tasman Front relaxed by the mid-Holocene, this subtropical isotopic influence on central Tasman Sea surface water δ 18 O could have diminished.
Taken together, a reduction in EAC/Tasman Front influences and increased STF influence by the mid-Holocene may have conjointly shifted central Tasman Sea δ 18 Osw signatures by about −0.2 to −0.3‰. That shift is close to the scale of change observed for Waitomo δ 18 O master values (early Holocene at 10-9 ka is −3.79‰; mid-Holocene at 6.5-5.5 ka is −4.07‰) when the speleothem data are Reduced EAC subtropical water inflow (warm, isotopically less negative) and weakening of Tasman Front influences between the early Holocene and mid-Holocene may be another key component to consider for Tasman Sea marine climate changes (and downwind connections to Waitomo). Modern ocean chemistry spatial relationships show a~0.2% northwest-southeast δ 18 O gradient for sea water situated between Australia and the North Island. Enhanced southward EAC transport during the early Holocene, with an expanded and invigorated Tasman Front, may have pushed isotopically less negative δ 18 O waters to the south that normally are positioned to the northwest of the North Island and adjacent to Australia's east coast. In this situation, isotopically less negative ocean waters could have been advected via westerly drift into the central Tasman Sea. When the EAC and inflow along the Tasman Front relaxed by the mid-Holocene, this subtropical isotopic influence on central Tasman Sea surface water δ 18 O could have diminished.
Taken together, a reduction in EAC/Tasman Front influences and increased STF influence by the mid-Holocene may have conjointly shifted central Tasman Sea δ 18 O sw signatures by about −0.2 to −0.3% . That shift is close to the scale of change observed for Waitomo δ 18 O master values (early Holocene at 10-9 ka is −3.79% ; mid-Holocene at 6.5-5.5 ka is −4.07% ) when the speleothem data are averaged at a similar resolution as the marine evidence. This hypothesis suggests ocean source effects may have contributed to Waitomo δ 18 O calcite changes between the early Holocene and mid-Holocene. An additional, intertwined component that may have helped to drive this type of ocean source region change is that the delivery of meteoric waters came from increased frequency of northerly quadrant airflow across lower latitude maritime ocean sources (particularly on a seasonal basis; discussed below).
Two of three available transient model counterparts (CCSM3 and ECHO-G; Figure 7, bottom row) that we examined cover the mid-to-late Holocene and show a winter Tasman Sea SST warming trend set within a wider southwest Pacific temperature rise. This trend suggests regional conditions during the mid-Holocene were locally cooler than pre-industrial times [120]. Regional downscaling of global climate models show mid-Holocene multi-centennial scale surface air temperatures for New Zealand were also at least 0.2 • C cooler in DJF through to JJA, due to circulation changes and reduced insolation [121]. SSTs were also 0.2-0.5 • C cooler in the Tasman Sea west of the North Island, and cooler temperatures existed over most of the North Island in the mid-Holocene [122].
The climate model-calibrated Waitomo δ 18 O ivc winter temperature warming trend from~5 ka to the present appears different from marine proxy evidence that suggests continuation of stable temperatures, or development of cooler temperatures up to the pre-industrial period (see Section 4.3). Putting any seasonal proxy sensitivities aside, the divergence between the marine and terrestrial data implies oceanic source region influences (as described above) on δ 18 O trends in the Waitomo isotope master record may have diminished relative to other drivers that took hold from the mid-Holocene onward.

Did Regional Ocean/Atmosphere Circulation Consistently Influence Waitomo δ 18 O Trends?
Early Holocene winter warmth and reduced subpolar influences around northern New Zealand that were succeeded by cooler conditions by the mid-Holocene may be explained by previous work that promoted "two circulation modes" [110]. An early Holocene mode has been characterized by increased penetration of subtropical flow and southward expansion of the STF that coincided with reduced subantarctic water influence on the South Island of New Zealand. A hypothesis for these early Holocene marine conditions (in addition to caveats about seasonal biases for marine climate reconstructions) is that they existed due to weaker southern hemisphere westerly winds (SWW), arising from a reduced latitudinal temperature gradient. Diminished westerly winds with an expanded southward position and increased strength of the subtropical high could have helped to push the EAC more vigorously at the same time as the STF was displaced further south. Evidence for a more southern position for the subtropical high is supported by well-stratified waters east of New Zealand, which may have been promoted by less vigorous mixing and lighter surface winds with increased passive penetration of warm subtropical waters bathing the northern North Island (see details in Section 4.3). The circulation pattern that followed and that was established by the mid-Holocene was characterized by progressive strengthening of the SWW with an increased influx of subantarctic waters, cooling south of the STF, and rebound of the southern margin of the STF due to an increased latitudinal gradient across the New Zealand sector [110].
Holocene insolation changes that influenced SSTs via enhanced southern mid-latitude surface meridional temperature gradients are proposed as a driver of a poleward shift and intensification of the SWW [123][124][125]. The tropospheric meridional temperature gradient and the surface meridional temperature gradient are also significantly correlated over the New Zealand sector [25]. This means upper atmosphere circulation shifts driven by latitudinal insolation variations can propagate through to the near-surface environment to produce surface temperature anomalies. Mid-Holocene orbital forcing of atmospheric circulation, via steepened tropospheric meridional temperature gradients, also include narrowing and intensification of the winter southern subtropical jet over parts of the southern lower middle latitudes with spatially variable SWW impacts (Figure 7). Transient model simulations indicate negligible change or a weakening trend for austral winter westerly flow across New Zealand from the mid-to-late Holocene (based on CCSM3, Echo G, and LOVECLIM; 3 of 4 transient models used in this study [120]). That evidence is congruent with regional-scale climate model analysis for New Zealand that evaluated winds, synoptic types, and climate regimes for the mid-Holocene relative to pre-industrial conditions [121,122] (see Figure A1 in Appendix C for synoptic types and climate regime details). Regional atmospheric circulation patterns across the hemisphere also show important longitude-and season-specific idiosyncrasies that arise from a wider SWW response to tropospheric temperature gradient changes [68].
A SWW strengthening and poleward shift can arise from a CO 2 increase, which has been demonstrated over the late 20th and early 21st century under a warming climate (e.g., [126,127]). Ice core evidence indicates increased atmospheric CO 2 from the mid-to-late Holocene occurred [128]. A mid-to-late Holocene central Tasman Sea climate state change may have therefore been partly driven by both greenhouse gas forcing and regional atmospheric circulation impacts on the surface climate.
Some of the associated regional changes in the transient model data also appears heterogeneous at larger spatial scales over the Southern Hemisphere through the Holocene, but these differences can be reconciled via understanding longitude-specific responses to SWW. We also note a subtle La Niña/positive SAM-like SST spatial pattern in models across the Pacific basin for the mid-to-late Holocene [120]. That type of climate pattern is typically associated with more frequent northerly and easterly quadrant flow with warmer regional temperatures across the North Island, which is consistent with tightening and strengthening of SWW to the south of New Zealand and an increase in subtropical inflow (at least on a seasonal basis). This pattern may suggest the increased importance of (and/or influence from) specific modes of variability that New Zealand is sensitive to (including ENSO [129], SAM [130] and Zonal Wave 3 [131]) that are known to work in conjunction [45,97]. However, further work with more finely resolved proxies would be required to investigate the potential role those climate modes had in guiding any observed changes during the Holocene.

Caveats for Interpreting δ 18 O and Reconciling Different Drivers of Holocene Changes
Spatially heterogeneous atmospheric circulation features across the Southern Hemisphere for the mid-Holocene, as seen climate model simulations, indicate reduced zonal flow across South American and Australian sectors of the high-middle and polar latitudes [68,132]. As a consequence, extratropical westerlies may have been weaker during the mid-Holocene over certain parts of the southern hemisphere, for certain seasons, despite a more intense overall upper tropospheric southern jet stream at times [68]. This example highlights different outcomes can arise from simultaneous SWW structural changes, and it emphasizes longitude-and latitude-dependent responses to prevailing atmospheric circulation changes (e.g., some places experience relative SWW weakening when other areas experienced relative SWW strengthening; see supporting data summarized in [110]).
Significant spatial differences that are observed for models and field evidence reinforces caveats from drawing broad conclusions about large scale processes (or about climate modes) from proxy records that capture local scale (or longitudinally constrained) conditions. This consideration reiterates commonly known limitations for proxy-model comparisons that can be surmounted by wider data networks. It also reemphasizes caution in drawing direct comparisons between Waitomo and other Southern Hemisphere sites further afield.
Nevertheless, Waitomo appears to be an important location within the wider Southern Hemisphere proxy network that records possible impacts from regional atmospheric circulation and Tasman Sea surface climate changes. Along with other New Zealand evidence, the new Waitomo δ 18 O ivc isotope master record can help to test hypothetical models of Quaternary climate change that consider how synoptic systems [21] and oceanographic boundaries are interconnected at the hemispheric scale via the SWW.

Conclusions
New Zealand speleothem data were age modeled using Bayesian principles with COPRA [31], and these new data are now held in the SISAL database for future use by the palaeoclimate research community. We briefly reviewed the national spatial and temporal coverage of isotopic series from New Zealand speleothem records now held in the SISAL database. The majority of our effort then focused on using COPRA and other techniques to develop a δ 18 O isotope master record for Waitomo that covered the Holocene, and interpreting that record within a proxy-model context. Several strengths and weaknesses exist for different master record construction approaches (see Section 4.2). From our findings, we suggest evaluating several master record methods in parallel while considering a range of factors (including the environments the speleothems grew in, the amount of local data replication, and the timespan speleothem collections cover) can help to determine whether the intended application of an isotope master record is appropriate.
Leading hypotheses for the millennial-scale trends in the Waitomo δ 18 O isotope master record suggests it reflects temperature (and isotopic) traits from different climatic sub-regions within and around the Tasman Sea. A key factor in our interpretation of this record is that the latitudinal temperature gradient and associated fluctuations in the SWW winds produce dynamic subpolar and subtropical impacts on regional atmospheric circulation, which in turn drive Tasman Sea and Waitomo surface climate conditions. During the early Holocene, a reduced latitudinal temperature gradient may have led to several impacts (including diminished SWW near New Zealand, displaced the STF southern boundary, increased Tasman Front influences, and enhanced EAC flow along the western Tasman Sea) that produced warmer winter temperatures, which are reflected in the model-calibrated Waitomo δ 18 O isotope master record.
The transition from the early-to-mid Holocene suggests the associated cooling of the southern high latitudes steepened the latitudinal gradient to the east of New Zealand [110]. Increased upper level and surface SWW brought a range of changes (including increased upper ocean mixing, increased northward penetration of SAF waters east of New Zealand, a northward rebound of the STF, and reduction of EAC inflow into the North Tasman Sea) that diminished subtropical water influences around northern New Zealand, and produced a cooling winter temperature trend for Waitomo.
During the mid-Holocene, we observe lowest δ 18 O ivc values, which we interpret as the coolest temperatures for JJA at Waitomo based on how we calibrated our record. Strong seasonality, and exclusion of frost tender taxa with increased frosts during winter (and flanking months) at this time are suggested from other New Zealand evidence [18]. Then, from the mid-to-late Holocene, CO 2 concentrations increased, which likely contributed toward a stronger and poleward-shifted SWW belt. Impacts from that atmospheric circulation change would have led to more mixing, coinciding with strong surface climate heterogeneities over portions of the mid-and-high latitudes and an increase in subtropical ridge presence over parts of the lower middle latitudes (25-40 • S).
The Waitomo Holocene δ 18 O isotope master record exhibits millennial-scale temperature trends with distinct periods of centennial-scale variability that can be corroborated using a range of available higher-resolution proxy evidence. However, we also recognize that the conceptual model of Holocene climate change offered in this study would greatly benefit from more detailed palaeoproxy data transects on the eastern (Pacific) and western (Tasman Sea) sides of New Zealand, in addition to a better understanding of the controls on precipitation δ 18 O variations over space and through time.

Future Work
New Zealand speleothems significantly contribute to palaeodata coverage for the Southern Hemisphere mid-latitudes. They offer a unique resource for defining a wide range of Late Quaternary natural climate and environmental variations. Regional richness of other proxy records in New Zealand also helps to support large-scale spatial interpretations from the karst record. Archived speleothems that are well-dated provide a valuable legacy resource for the palaeoclimatology community that can be re-assessed with new analyses (e.g., trace element and physical properties) to complement interpretations from extant stable isotope records. Potential to build on the extant stable isotope data foundation requires ongoing stewardship of speleothem samples and the ability for researchers to access them.
Based on the New Zealand data lodged with SISAL so far, and what is forthcoming, it is clear many caves can be revisited in order to augment sample depth, apply new and more precise dating approaches, and improve climate-proxy calibrations. Enhancing our understanding of New Zealand speleothems that cover the last glacial cycle and beyond also requires high-resolution, long-term cave monitoring to support the development of quantitative climate transfer functions. Application of recent advances that have applied synchrotron-derived layer counting [133] paired with high-precision radiometric dating appear to be critical for this effort. Re-expansion of cave climate monitoring experiments where speleothems can be collected will help to evaluate how inter-annual to decadal-scale climate variability and extreme events (droughts, pluvials, fires) may be imprinted during carbonate deposition.
Northland (northern North Island), Hawke's Bay (eastern North Island), and north Canterbury (eastern South Island) are high priority for development of new speleothem data to fill current spatial gaps. The area north of Auckland including the cave systems near Waipu (both in northern North Island) have potential to provide speleothem material for palaeoclimatology. Notably, this region also has many palaeoproxies, including annually-resolved sub-fossil kauri tree rings that cover parts of the Holocene, the late glacial, most of MIS3, and also parts of MIS5 [134]. However, LGM depositional hiatuses or unconformities in some sedimentary records there are noted [135]. Thus, Northland speleothems may be able to provide continuous evidence of environmental changes spanning the last glacial-interglacial cycle relative to what has currently been reported from that region. Several Hawke's Bay speleothems are already archived, and additional replication from that area is warranted because of regional impacts from westerly circulation [21,40] variability and apparent drought sensitivity. For north Canterbury, which is also a relatively dry region in the lee of the Southern Alps, the situation is similar; but the challenge may be to find suitable speleothems (especially those that are in isotopic equilibrium with the local environment).
Application of COPRA for consistent age modelling aided our comparison of isotope master record techniques. Success or failure for splicing speleothem time series into a long master record is likely related to the amount of data available at a sub-regional level. This assertion is supported by our efforts to combine isotopic data starting within a cavern and working outward to additional data from neighboring caves within the same karst block. While data aggregation processes used to create isotope master records could introduce data artefacts when combining or adjusting series, increased replication has benefits for evaluating the timing and scale of relative trends and changes through time. Additional proxies (e.g., trace elements) derived from well-dated speleothems that are aligned on a common time scale would provide more temporally-specific details that could enhance isotope master record construction practices, in addition to aiding the re-assessment of individual speleothem sensitivities to climate (and in particular, hydroclimate).
Future upgrades to the wider Australasian speleothem data network are expected to enhance understanding of heterogeneous patterns of change driven by global climate teleconnections. Demonstrable links for many Australasian caves to convergence zone activity, temperature, and rainfall processes means they remain important pre-instrumental resources for interrogating planetary ocean-atmosphere circulation dynamics. Given the limited instrumental data coverage for much of Australasia, including New Zealand, it is clear speleothem perspectives can provide highly valuable, idiosyncratic information that addresses spatiotemporal biases or gaps for past environmental histories that are not covered in the literature. Funding: This work was undertaken as part of NIWA Strategic Science Investment Fund projects "Regional Climate Modelling" and "Climate Present and Past" which are funded by New Zealand Ministry for Business Innovation and Employment.

Acknowledgments: SISAL (Speleothem Isotope Synthesis and Analysis) is a Working Group of PAGES (Past Global
Changes) program. We thank PAGES for their support for this activity. We also thank the New Zealand Department of Conservation and all landowners for access to caves and permissions for sampling speleothems. Andrea Tuohy and Tom Whittaker are thanked for access to speleothem data. Bedartha Goswami is thanked for discussions about COPRA and conservation of age model output that assisted our aims, and Sebastian Breitenbach is thanked for training and assistance with COPRA. Duncan Ackerley, Sebastian Wagner and Johann Jungclaus are thanked for discussions and generously providing access to climate model simulation data. We thank the World Karst Aquifer Mapping project (WOKAM) team for data underpinning the karst region map in the Figure 1 inset. We are grateful for in-depth comments from Emily Judd, Pauline Treble and Kale Sniderman which greatly improved this manuscript. Laia Comas-Bru is especially thanked for supplying the inset in Figure 1 and for providing constructive comments on much earlier versions of this manuscript, in addition to her extensive coordination efforts related to SISAL.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A. Physical Geography of Karst Terrain of New Zealand
Regions covered in this work have been previously described by some of the authors in a technical report [47]. We reiterate that material for the benefit of researchers who wish to understand more about the local geology, climates, and past speleological investigations.
The Waitomo district is located at latitude 38.3 • S about 35 km inland from the west coast of the central North Island. Karst in this district (sometimes referred to as King Country karst) is principally in Orahiri Limestone of Oligocene age and Otorohanga Limestone of Oligocene-Miocene age, and rests unconformably on an impervious basement that is overlain by relatively silty and sandy clastic beds [24,37]. The karstified limestone is discontinuous and covers an estimated 800 km 2 in subcrop and outcrop. The area has a temperate oceanic climate with a mean annual temperature of about 13 • C and an annual rainfall averaging 1618 mm at Waitomo, although this increases to over 2300 mm some 10 km to the west. The average annual temperature of caves in the Waitomo area is about 12.8 • C. Podocarp-hardwood evergreen forest formerly covered the district [136], but much has been cleared for agricultural grassland in the last 100 years. Further background details and information on the regional context of the site are provided in previous work [12]. Speleogenesis is understood to have started mainly in the Pleistocene. Hence the upper levels of the longest cave in the region (Gardner's Gut) have speleothems dating from the Last Interglacial. However, it is likely that the formation of most caves in the area has taken place within the last half million years and often much more recently although some caves are >1.0 Ma. The six Waitomo speleothem chronologies used in this synthesis report for the Holocene analysis include two from Gardner's Gut Cave (GG1, GG2), three from Ruakuri cave (RK-A, RK-B, RK-C) and one from Max's Cave (Max's). Gardner's Gut Cave has a total passage length of 12.2 km and is located 3.5 km west of Waitomo village beneath an undulating karst surface varying from 120 to 180 m above sea level. Max's Cave is located about 8 km west of Waitomo. The Gardner's Gut speleothems were all obtained from approximately the same elevation (ca. 100 m above sea level), whereas the speleothem from Max's Cave was obtained somewhat higher (ca. 325 m above sea level).

Appendix A.2. Whakapunake Karst, Hawkes Bay
The Whakapunake area is approximately 29 km north-east of Wairoa in the Hawkes Bay District. Limestone landforms are well distributed along this eastern side of the country and are mainly of Pliocene to early Pleistocene age [37]. Typically, these limestones are young and poorly cemented, so caves and sampling opportunities are relatively limited, though these young porous limestones still contain caves. However, the limestone at Whakapunake is dense and crystalline [137]. Streams have intensely dissected the coverbeds and soils over the limestones. Overlying vegetation is mixed podocarp, broadleaf and hardwood forest and includes many tall trees species [138], although on Whakapunake most natural forest has been cleared and has been replaced by pine plantations. Steep terrain in this region has developed due to regional tectonic uplift along the western side of the Hikurangi-Tongan trench subduction zone. Annual rainfall averages some 1550 mm, with average annual temperatures about 12.5 • C. The speleothem chronologies from the district used in this synthesis report comprise speleothems from Disbelief and Te Reinga caves. Both these caves are located in the Wairoa River catchment near the town of Te Reinga [139]. Disbelief Cave is at an elevation of about 600 m and is approximately 120 m long. Te Reinga Cave is 2690 m long and 134 m deep (labelled Te Ranga in [139]). Speleothems were sampled in Te Reinga at about 100 m above sea level (asl).
Westland karst [140]. The limestones tend to be of variable thickness and lithology-ranging from thin and sandy to crystalline [140]. Between Charleston and Punakaiki the karst extends inland for several km's uplifted in places to 430 m asl. The West Coast has a mild maritime climate with offshore waters being transitional between subtropical and sub-Antarctic. Sea surface temperatures range between a mean winter minimum of 13.3 • C and a mean summer maximum of 18.4 • C. Punakaiki and Westport have mean annual temperatures of 13.7 • C and 12.3 • C respectively, but temperatures in nearby caves have been measured at 10-10.1 • C. Mean annual precipitation is 2346 mm at Punakaiki, but increases inland with altitude on seaward slopes. Most precipitation falls as rain, but at sites above 400 m precipitation can fall as snow in winter. Prevailing onshore westerly winds convey a strong moderating effect on coastal air temperatures. The area is covered with rainforest dominated by podocarp and evergreen hardwood-broadleaf species. Three stalagmites near Punakaiki were used in this synthesis report-one from Babylon Cave (BN1), one from Hollywood Cave (HW1) and another from Wazpretti Cave (WP1). Hollywood Cave is 4.3 km long. It is located near Charleston on the West Coast of the South Island, between Westport and Punakaiki. The entrance to the cave is 5 km inland from the coast and is located at the bottom of a doline (~130 m) set in an undulating plateau with summits to 380 m above sea level. HW1 speleothem is a stalagmite of simple elongate morphology and is just over 500 mm long. It was situated in a link passage between Streamway One and Streamway Two about 200 m into the cave beyond the restricted entrance passage and at least 50 m below the undulating epikarst surface. Babylon Cave is located near to Fox River and is approximately 800 m long. Wazpretti Cave is located about 6 km inland along Bullock Creek and at 90 m asl. Speleothem samples were obtained in Babylon Cave at 175 m and in Wazpretti at 95 m asl.
Appendix A.6. Te Anau, Fiordland Karst in the Te Anau area is developed in Tunnel Burn Formation, an Oligocene bioclastic limestone 30-80 m thick [37]. The Aurora Cave system is located along the western side of Lake Te Anau, the largest lake in South Island. It comprises an upstream section called Aurora Cave separated by a sump from a shorter downstream resurgence section named Te Ana-au Cave (a tourist cave). Aurora Cave was formed by the stream draining Lake Orbell, in a hanging tributary of Lake Te Anau [11]. The Lake Orbell outflow stream (Tunnel Burn) descends 267 m through the cave to Lake Te Anau [11]. Aurora Cave has 6.4 km of passages. Speleothem deposits in Aurora Cave are particularly important because they have shed light on the timing of glacial advances and interstadial periods [11]. Aurora 3 speleothem was obtained from about 260 m above sea level in a downstream part of the cave. Calcite Cave is located some 13 km south-west of Aurora Cave at 1025 m elevation on Mt. Luxmore. It is above the tree-line and beneath alpine herb-field. The cave ascends from its entrance and samples were obtained at about 1035 m. The cave is within the Tunnel Burn Formation and approximately 200 m long. Climate of this south-western area of the country is temperate and dominated by the prevailing westerlies. Air masses from the tropics and polar regions can reach this region, with heavy rainfalls and cold showery conditions not uncommon. At an elevation of 215 m asl the Te Anau climate station has a mean annual temperature of 9.3 • C, although clearly elevation, exposure and proximity to the sea will affect the temperature of particular localities across this area. Mean annual precipitation is 1180 mm at Te Anau, with a sharp westward increase in rainfall.
Appendix A.7. Doubtful Sound, Fiordland Doubtful Xanadu Cave is located above Kellard Point in Doubtful Sound, Fiordland to the northwest of Lake Chamberlain. The stream-sink cave entrance (960 m above sea level) is situated immediately above tree line, with overlying vegetation mainly comprised of alpine tussock. Glacial geomorphic features in close proximity to the cave clearly indicate a landscape history deeply intertwined with extensive glacial activity. Local bedrock is composed of a salmon pink marble that is interbedded with dark grey to black schist, which gives an unusual horizontal striped appearance to the country rock that comprises this cave. Meta-sedimentary strata dip toward Kellard Point. There is approximately 700-800 m of passage in Doubtful Xanadu, but the main passage distance is 360 m long with a vertical descent of 50 m. Regional tectonic activity has fractured the country rock into multiple sets of joints that have been exploited by groundwater to form passages with an average of 1-3 m in width for this cave system. The speleothem sample was obtained from a passage at about 950 m altitude. There are no climate stations in Doubtful Sound, and the closest meteorological record is from Milford Sound, which extends back to 1927 AD. There are significant altitude and latitude gradients in precipitation due to mountainous terrain in this region. The New Zealand virtual climate station network (which estimates climate conditions based on interpolations within a network of stations) suggests the Doubtful Sound long-term median rainfall exceeds 4 m per year, with more precipitation during summer and autumn than winter and spring. Median annual temperature is close to 7 • C, with a summer median of 11 • C and a winter median of~4 • C. The area receives approximately 1450 bright sunshine hours per year, and has a median annual wind speed of~6 m/s [25]. Influences on New Zealand speleothem calcite δ 18 O composition have been partitioned into drip water (seepage water) effects and cave temperature effects (related to thermodynamic fractionation) [25]. In the context of New Zealand speleothem environmental monitoring [60], few observations are published on infiltration, precipitation residence time in the epikarst, karst aquifer storage processes, mixing of groundwater and subterranean flow styles (diffuse flow vs. fracture flow; e.g., [16]), but a recent contribution has significantly expanded our understanding of these factors in the Waitomo region [34]. All of those local aspects can potentially contribute toward enrichment or depletion of 18 O of seepage water and therefore speleothem calcite δ 18 O composition. It is likely that some of those processes also cause variations between speleothem δ 18 O within the same cavern and between caves in the same region.
Rainfall amount effects on δ 18 O of precipitation (δ 18 O p ) shows a seasonal effect exists across the year [12], with less negative δ 18 O p during summer when rainfall is reduced and more negative δ 18 O p when rainfall is higher during winter [60]. Precipitation amount effects on δ 18 O is expected to be larger for warmer climates than cooler temperate climates, and all data in this synthesis are drawn from sites with mean annual temperatures below 15 • C [60], though evidence from the Waitomo region suggests an amount effect in seasonal precipitation δ 18 O variations [34]. In addition, rainfall amount effects on δ 18 O p from distinct storms (with increased negative values arising from Raleigh distillation) are expected to be smoothed out in low-frequency speleothem analyses (>10 y), which are typical resolutions for New Zealand speleothem time series (see Table 1).
The Tasman Sea region (upwind of prevailing flow that blankets both main islands) is also expected to influence local δ 18 O signatures in New Zealand speleothems [47]. Within New Zealand, terrestrial temperature regimes reflect the traits of proximal oceanic waters that air passes over for 1000 s of kilometers before it arrives at a site [144]. On glacial-interglacial time scales, there were latitudinal shifts of the Subtropical and Subantarctic Fronts (STF, SAF) that impacted on regional sea surface water traits ( Figure 1). The location of oceanic fronts west of New Zealand can impart a significant influence on terrestrial air temperature patterns [82,110]. Sea surface temperature and sea water salinity attributes of the subtropical and subpolar water sources that characterise the STF and SAF are intertwined in δ 18 O sw , which influences δ 18 O p . Thus, palaeoceanographic studies are useful complementary data that can assist interpretations of δ 18 O in New Zealand speleothems. Because of the relatively slow inertia of oceanic changes, we recognise that important glacial-interglacial changes for New Zealand speleothem δ 18 O may have apparent regional differences due to oceanic circulation (and marine front boundary location changes). It is also for these reasons that pushing a Holocene-based transfer function in to the late glacial and the LGM has not been attempted in this study.
Precipitation δ 18 O analyses [145] show increasingly more negative δ 18 O from northern to southern New Zealand based on annual average δ 18 O p . Analyses on modern straw calcite from caves reflects this latitudinal gradient [12,22]. Modelling of precipitation isotopes on a global scale [146] also shows New Zealand has a slightly larger intra-annual range (termed the seasonal effect [12]) for δ 18 O p in the north (annual mean of about 5.5% at Kaitaia; −2.3% in January to −7.8% in July) relative to the south (annual mean of about 4.5% at Invercargill; −5.1% in February to −9.7% in July). Based on the mean annual temperatures for both these sites, there is a positive relationship between δ 18 O p and temperature on the order of 0.47% / • C [12], but these relationships are slightly different through the year for Tmax, Tmin and Tmean. A check of that assertion using Online Isotopes in Precipitation Calculator [146] closely agrees (0.41% / • C). Monthly isotopic information based on OIPC [146] and local Waitomo temperature data, however, show a relationship of 0.54% / • C which is situated more inland than Kaitaia or Invercargill. This relationship has a shallower slope for winter (June-August 0.32% / • C) when data from Kaitaia, Waitomo and Invercargill are examined in combination, but overall reiterates the point that δ 18 O p and temperature for New Zealand are positively correlated.
Previous work sampling rainwater and dripwater at Waitomo has shown δ 18 O p intra-seasonal isotopic ranges can be large (>6% ) but variability for temporally-aligned cave dripwater δ 18 O samples is very small (<0.5% ) [34,60]. This finding suggests cave seepage waters at Waitomo are thoroughly mixed before they reach the cave drip point, and also that the drip water source above the cave is potentially of sufficient size to attenuate intra-seasonal variability. Therefore, there should be a minimal influence for precipitation seasonality on low frequency signals that are the focus of the Waitomo speleothems examined in this study, and the δ 18 O of drip water at that site likely represents the mean annual precipitation value. We refer readers to the most recent work for additional details [34].
Along with different precipitation source regions, the regional temperature anomalies that arise from inter-seasonal synoptic weather variability [21,40,41,45] can impact on cave climates. The δ 18 O measurements from speleothems that are taken to represent average climatological conditions (with resolution dictated by sampling intensity, speleothem deposition rate and chronology constraints) may also contain an aspect about the persistence of regional climate regimes that convolve specific source region δ 18 O sw and also the temperature conditions of those source areas.
Speleothem δ 18 O variability that arises from the regional atmospheric circulation may dominate the second-order oxygen isotope trends at multi-decadal to centennial scale [21]. However, it has been noted that there are competing, opposite effects on speleothem δ 18 O from a range of sources [12,147] that can mute or reinforce the isotopic signal related to ambient temperature fluctuations. In the New Zealand situation, the cave temperature effect has been shown to influence the δ 18 O signal in speleothem calcite [22]. Previous work has indicated opposite signs exist for the cave temperature effect (negative) and the relationship between ambient temperatures and δ 18 O p (positive). The equilibrium fractionation factor (the cave temperature effect) for speleothems examined in this study is expected to range between −0.23% / • C near Waitomo and −0.26% / • C in Fiordland based on the mean annual temperature of each site. In a maritime climate like New Zealand, depending on what basis one starts from for the isotopic composition of seasonal recharge waters (which may strongly influence the δ 18 O signal following recent assertions [65]) and the equilibrium fractionation factor, the directionality of past climate change from δ 18 O will be similar and should be positive (i.e., less negative δ 18 O is likely to be warmer, more negative cooler), but the absolute scale of temperature reconstructions may be different. Tests of temperature reconstruction efficacy with and without the equilibrium fractionation factor used and held within a multiproxy framework may therefore be warranted.
When the counteracting effects of equilibrium fractionation (ca. −0.23% / • C) and temperature (0.4 to +0.54% / • C) on precipitation are summed (+0.17 to 0.31% / • C), a range of +/−1.5% in speleothem calcite δ 18 O relative to millennial-scale base climate state average conditions is probably realistic to attribute to temperature variability driven by regional atmospheric circulation operating on decadal to multi-centennial scales.

Appendix B.2. δ 13 C in New Zealand Speleothems
Quantitative transfer functions between δ 13 C and meteorological data are not numerous or based on long observation series [60], so many relationships for New Zealand speleothems are theoretical. It has been postulated that soil CO 2 influences speleothem δ 13 C variability via C 3 vegetation and forest productivity changes [30]. This assertion follows previous work that suggests diminished water flow rates enrich δ 13 C, and explains why more negative Holocene δ 13 C values that contrast with more positive LGM values (~8% difference) arise from interglacial-glacial differences in dryness and vegetation cover [87,148].
Nettlebed cave speleothem trace element chemistry, ultraviolet luminescence intensity (UVL), and Uranium isotopic data support the interpretation of δ 13 C as a palaeohydrology proxy. Millennial-scale UVL variations closely mimic coeval δ 13 C patterns and have been hypothetically linked to seepage water organic acid content flux that impacts speleothem fluid inclusion density, an assertion now supported by experimental evidence [149]. In this situation, increased (or decreased) speleothem growth rates that contribute to calcite formation should align to increases (or decreases) in percolation water flow rates, elevated (or reduced) biological activity in and above the epikarst, and warmer (or cooler) temperatures. It has also been outlined by previous workers [52] how increased groundwater residence driven by reduced effective precipitation above the cave drives elevated magnesium concentrations in speleothems [150]. From this result, and because reduced changes in magnesium and other trace elements closely parallel δ 13 C patterns, both types of records are interpreted as driven by effective precipitation changes on glacial-interglacial timescales [150].
Williams et al. [22] concurred that δ 13 C variations capture the combined effects of biological activity and local water balance changes. More positive δ 13 C values that coincide with the global LGM are associated with dry conditions, while more negative values after the last termination are associated with elevated biological activity and wetter conditions [22]. Interpretation of a network of speleothem δ 13 C signals across New Zealand's heterogeneous climatic regions using an atmospheric circulation regime paradigm also suggests the water balance link of δ 13 C to effective precipitation is plausible [25,40]. In periods when precipitation is high, the flushing rate in the epikarst has higher potential to increase in frequency, and therefore the residence time of water above the cave is expected to decrease [24]. This decreases the potential for inorganic δ 13 C enrichment of δ 13 C calcite , suggesting wetter periods should have more negative δ 13 C values. With this in mind, this interpretive framework could be validated using a combined δ 13 C, 14 C approach [151].
Water balance in the epikarst, atmospheric CO 2 , and vegetation composition changes over the cave (C 3 /C 4 pathways; see comments in Figure 3) are important to also consider for interpreting New Zealand speleothem δ 13 C. This type of proxy record may also have more important links to local or global temperature change in certain situations. Some δ 13 C signatures may not reflect only one aspect of climate strongly, or be affiliated with only one environmental element (like vegetation type), so interpretations within a multi-proxy framework are recommended [21]. It is clear that determining the relative importance of global and local factors (CO 2 , local vegetation type changes, whether a system is open or closed, climate conditions, epikarst processes like flushing) at any site are required to evaluate what δ 13 C controls are most prevalent.

Appendix C. Details on Regional Synoptic Weather Types and Weather/Climate Regimes of New Zealand
Quaternary 2020, 3, x 31 of 39 Appendix C. Details on Regional Synoptic Weather Types and Weather/Climate Regimes of New Zealand Figure A1. (Top) New Zealand synoptic weather types related to three key regimes (Trough, Zonal, Blocking; figure modified from [21,25,40,41,50]). Synoptic weather types associated with each regime are as follows: Trough group-T, SW, TNW, TSW; Zonal Group-H, HNW, W; Blocking group-HSE, HE, NE, HW, R (percentages indicate climatological occurrences in a twice daily k-means clustering of MSLP using the NCEP1 reanalysis. (middle) Example spatial patterns for three weather/climate regimes used to group the New Zealand synoptic types via similar rainfall and temperature outcomes (bottom). Details for New Zealand average associated spatial conditions arising from the synoptic types in each of the regimes for rainfall and temperature, and association