Evaluation of PMIP2 and PMIP3 simulations of mid-Holocene climate in the Indo-Pacific, Australasian and Southern Ocean regions

Abstract. This study uses the simplified patterns of temperature and effective precipitation approach from the Australian component of the international palaeoclimate synthesis effort (INTegration of Ice core, MArine and TErrestrial records – OZ-INTIMATE) to compare atmosphere–ocean general circulation model (AOGCM) simulations and proxy reconstructions. The approach is used in order to identify important properties (e.g. circulation and precipitation) of past climatic states from the models and proxies, which is a primary objective of the Southern Hemisphere Assessment of PalaeoEnvironment (SHAPE) initiative. The AOGCM data are taken from the Paleoclimate Modelling Intercomparison Project (PMIP) mid-Holocene (ca. 6000 years before present, 6 ka) and pre-industrial control (ca. 1750 CE, 0 ka) experiments. The synthesis presented here shows that the models and proxies agree on the differences in climate state for 6 ka relative to 0 ka, when they are insolation driven. The largest uncertainty between the models and the proxies occurs over the Indo-Pacific Warm Pool (IPWP). The analysis shows that the lower temperatures in the Pacific at around 6 ka in the models may be the result of an enhancement of an existing systematic error. It is therefore difficult to decipher which one of the proxies and/or the models is correct. This study also shows that a reduction in the Equator-to-pole temperature difference in the Southern Hemisphere causes the mid-latitude westerly wind strength to reduce in the models; however, the simulated rainfall actually increases over the southern temperate zone of Australia as a result of higher convective precipitation. Such a mechanism (increased convection) may be useful for resolving disparities between different regional proxy records and model simulations. Finally, after assessing the available datasets (model and proxy), opportunities for better model–proxy integrated research are discussed.


Introduction
Proxies give indications of past climatic conditions 1 and can be used to assess the ability of atmosphere-ocean general circulation models (AOGCMs) to represent past climate states. Moreover, if past climate states can be reproduced adequately with AOGCMs, then there can be a degree of confidence in their ability to simulate future climate change. A direct comparison between model and proxy data is difficult, particularly when AOGCM grid spacing is ∼ 100 km (or more) and proxies represent climatic information at a specific place or occasionally over a broader region. A method of bridging this scale gap is to "upscale" various regionally coherent proxy reconstructions to a spatial scale that is resolved by AOGCMs. Conversely, large-scale AOGCM data can be "downscaled" using known circulation characteristics to provide local-scale (≤ 100 km) estimates of climatic variables (e.g. precipitation and temperature). Such an upscaling approach was adopted in Lorrey et al. (2007Lorrey et al. ( , 2008 to use regionally coherent climate proxy data (temperature and precipitation) to infer circulation characteristics over New Zealand. Subsequently, Ackerley et al. (2011) applied the reverse method to downscale coarse-resolution AOGCM data to infer regional temperature and precipitation characteristics over New Zealand, which provided a platform to evaluate the merits of both model and proxy datasets. The opposing approaches, i.e. upscaling the proxies versus downscaling the models, employed by Ackerley et al. (2011) and Lorrey et al. (2007Lorrey et al. ( , 2008, enabled both datasets to be compared in a more direct way and provided a platform from which to investigate their respective merits and shortcomings. Such upscaling/downscaling approaches for comparing proxy and model data in a direct and meaningful way have not been attempted yet over the broader Australasian sector. The southern Maritime Continent, Australasia and Southern Ocean (immediately due south of Australia) regions considered by the Australian INTegration of Ice core, MArine and Terrestrial records (from now -OZ-INTIMATE) span from 10 • N to 60 • S and 115 to 155 • E (Fig. 1). The regions incorporate the tropical, arid, temperate and southern Indian Ocean/Southern Ocean climatic zones (see Fig. 1). While climate reconstructions from proxy records for the last 35 000 years are discussed individually elsewhere (see Bostock et al., 2013;Petherick et al., 2013;Reeves et al., 2013a, b, for more details), it is the schematic representation of "simplified patterns of temper-ature and effective precipitation 2 between regions" (p. 24 of Reeves et al., 2013a, their Fig. 4) that provides the opportunity for direct model-proxy comparison in an up-scaled approach. In particular, this approach allows us to move away from the purely descriptive and towards testing some of the dynamical mechanisms responsible for the proxy response, providing a deeper understanding of some of the key drivers of change. While reference is made to OZ-INTIMATE (as many relevant papers were produced through that initiative) this research is a contribution to the Southern Hemisphere Assessment of PalaeoEnvironment (SHAPE) initiative, an INQUA International Focus Group (project 1067P). SHAPE is the successor to the Australasian INTIMATE project, which broadens the geographical scope to encompass the entire Southern Hemisphere. One of the key remits of SHAPE is to develop model-proxy comparisons to both help understand the dynamical mechanisms behind past changes and also to test the robustness of the models. This paper presents such a comparison, which may provide a starting point for further model-proxy studies within the SHAPE initiative.
Regional reconstructions of temperature over the southern Maritime Continent, Australasia and the Southern Ocean (immediately due south of Australia) have previously been compared with Paleoclimate Modelling Intercomparison Project (PMIP; Joussaume and Taylor, 2000) AOGCM simulations for the past millennium (PAGES 2k-PMIP3 group, 2015). The PAGES 2k-PMIP3 group (2015) study provides an estimate of mean temperatures for the whole region (specifically 0-50 • S and 110-180 • E) instead of subdividing into smaller climatic zones (e.g. tropical versus mid-latitude or continental versus maritime). Here, the focus is on the mid-Holocene (ca. 6000 years before present; 6 ka) experiment of PMIP2 and PMIP3 (see Braconnot et al., 2007Braconnot et al., , 2012, in comparison with the pre-industrial era (ca. 1750 CE). Other studies have also evaluated data from the PMIP 6 ka simulations regionally, for example over South America (Prado et al., 2013) and Europe (Brewer et al., 2007;Mauri et al., 2014). This study therefore provides the first instance (to our knowledge) of this being done for the 6 ka time slice, which focusses on the Australasian region. The 6 ka time slice is a fairly "unremarkable" period in the climatic history of Australasia, but as a result there is an opportunity to evaluate the impact of weaker seasonality (lower austral summer insolation) on the climate. Furthermore, this paper presents a case for the development of more seasonally resolved proxy reconstructions to compare with the models. This study also provides an overview of the climatic conditions in the Australasian region at 6 ka from the PMIP archive, which was not discussed in the Braconnot et al. (2012) summary.
The aims of this study are fourfold. First, the study provides a rigorous assessment of the temperature, precipita-  Figure 1. A map of the geographical region under consideration in this paper. Overlaid are the borders and names of the regions referred to in the text and correspond to those of Reeves et al. (2013a). The divisions are also broadly consistent with the climate classifications (Köppen-Geiger) described in Peel et al. (2007). tion and circulation characteristics of the PMIP 6 ka experiments (Braconnot et al., 2007) over the regional domains outlined by Reeves et al. (2013a) (see Fig. 1), where proxy data display some spatial coherence. Second, the paper shows where the models and proxies (adapted from Reeves et al., 2013a, and references therein) agree or disagree by comparing them directly against each other. Third, where the models and proxies agree, the processes that are responsible for the climatic state are highlighted. Conversely, where the models and proxies conflict to some degree, the fourth aim is to provide an explanation as to why any dispute (or uncertainty) arises. It is not the intention of this work to say the models or the proxy interpretations are incorrect; instead the intent is to show that proxy-model agreement gives confidence in our assessment of past climate and the dynamical mechanisms behind the changes recorded by the proxies. Conversely, disagreement provides a key opportunity to re-focus our efforts and resolve the issue in an integrated way.
An overview of the data and methods used in this analysis are presented in Sect. 2. A synthesis of the model-simulated and proxy-inferred climates for the regions in Fig. 1 at 6 ka are presented in Sect. 3, which also includes (i) a discussion of the processes responsible for the climatic state in each region when the models and proxies agree and (ii) a consideration of the limitations of the models/proxies where there is some inconsistency or uncertainty.
Suggestions as to where future efforts should be focussed are presented in Sect. 4 and a summary of the main results and conclusions are given in Sect. 5.

Proxy data
Much of the palaeoclimate data used to support the analysis in this paper are derived from regional reviews of Australasian past climate Petherick et al., 2013;Reeves et al., 2013a, b), which had the initial purpose of developing a climate event stratigraphy for that region as a contribution to the previous INTIMATE program. The selection criteria for INTIMATE records were as follows: i. that they are continuous and cover the period of interest; ii. that they retain centennial to millennial scale resolution; iii. that they have robust chronologies.
As this is rarely achieved in Australia, the OZ-INTIMATE synthesis included discontinuous records that are well dated (given the limitations of the available methodologies and datable material) and include key intervals of change with a reconcilable proxy response to climate. A combination of high-resolution, centennial or better scale reconstructions (e.g. marine, speleothem, coral records) is included along with discontinuous geomorphic records (e.g. fluvial, lake shore, dune, glacier) and well-constrained qualitative and semi-quantitative biological records (pollen, diatom, ostracod, charcoal and geochemistry). Note that records that have been published since 2013 that fit these criteria have been included in our present study, where they meet the necessary criteria and add clarity to previous interpretations. Whilst the original interpretations of the records were maintained, the context of the site, limitations of the proxies and chronological integrity of the records were also considered. A database of the records included in this paper can be found at https://doi.org/10.1594/PANGAEA.879515. Included in the database are details of the original records upon which this and the OZ-INTIMATE compilations Petherick et al., 2013;Reeves et al., 2013a, b) are based. Where possible, location, proxy, climate sensitivity and dating and sampling resolution are included, as well as links to source data storage. A shorter document, containing all the proxy data that are referenced directly in this article, is also provided as a Supplement for this paper. Reeves et al. (2013a) subdivide the Australasian-southern Maritime Continent region into four main zones, which are indicated by the solid lines in Fig. 1, namely tropical, subtropical, temperate and Southern Ocean. Reeves et al. (2013a) also subdivide those four regions into the perennially wet equatorial tropical north-west and north-east (termed TNW and TNE, respectively); the monsoonal tropical southwest and south-east (termed TSW and TSE, respectively); the arid subtropics (termed StA) in the continental interior of Australia eastward of the Great Dividing Range; the temperate east and south (termed TeE and TeS, respectively) adjacent to the coast with maritime climates; and the northern and southern Southern Ocean (termed NSO and SSO, respectively). These regions broadly correspond with the Köppen-Geiger climate classification system (Peel et al., 2007) and correspond with regions where groups of proxies display their strongest agreement spatially. Changes in temperature and precipitation for different time slices were presented in Reeves et al. (2013a) relative to the previous time slice for each of these subregions. For example, if a subregion is warmer (cooler) in one time slice relative to the previous, the whole box is shaded red (blue). A similar analysis and colour scheme is also applied to effective precipitation. Such an analysis can be applied easily to surface temperature and precipitation data from AOGCMs by area averaging over the regions in Fig. 1. The area-averaged temperature and precipitation fields between different periods could then be quantified in the models and compared directly with the proxy data. The OZ-INTIMATE compilations focussed on trends from a previous period to the next and generally grouped the interval around 6 ka as part of a broader mid-Holocene phase or, distinct from the 8-7 ka period with respect to temperature and precipitation conditions Reeves et al., 2013a). Therefore, in this study, the reconstructions assessed in Reeves et al. (2013a) (and references therein) are re-evaluated to provide an estimate of the climatic state around 6 ka (±500 years) relative to the pre-industrial era in order to make a more direct comparison with the AOGCMs (Sect. 2.2).

Model simulations and boundary conditions
There is a vast amount of palaeoclimate AOGCM output that is freely available from PMIP (Joussaume and Taylor, 2000;Braconnot et al., 2007Braconnot et al., , 2012. The PMIP initiative includes data from transient simulations of the last millennium along with time slice simulations of the pre-industrial era ca. 1750 CE (0 ka), the mid-Holocene (6 ka) and the Last Glacial Maximum (21 ka). In this study we make use of coupled AOGCM simulations conducted for Phases 2 and 3 of PMIP (PMIP2 and PMIP3, respectively) for 6 ka. Full details of the experiments run, and evaluations of the simulated responses can be found in Braconnot et al. (2007Braconnot et al. ( , 2012, Harrison et al. (2014) and Taylor et al. (2012).
Output from the mid-Holocene (6 ka) and the preindustrial control experiments (0 ka) is used here. The boundary conditions (for example, orbital parameters and greenhouse gas concentrations) are given for the 6 and 0 ka simulations in Table 1. The available simulation data from PMIP2 (both 0 and 6 ka) and PMIP3 (6 ka) are 100 years in length and all years are used in the analysis below. The 0 ka simulations from PMIP3 vary in length from 100 to 1000 years, but all available years are used in the analysis below to minimise the impact of model internal variability. In all, 32 different model simulations (18 from PMIP2 and 14 from PMIP3) are used in this study, which are listed in the Supplement (Table S1), as individual models are not considered in this study. The original grid configurations of the models (along with the relevant references for each model) are given in Table S1; however, all model data were bilinearly interpolated to a common longitude-latitude grid (2.5 • × 2.5 • ) before undertaking the analysis below for ease of comparison. Moreover, due to the large number of models considered in this study, a detailed analysis of the individual model performances is not undertaken. The specific details of the individual model performances are discussed in the relevant references given in Table S1.
In order to show the impact of the orbital parameter changes, the zonal-mean change in incoming solar radiation at the top of the atmosphere (insolation) is plotted in Fig. 2a for 6 ka relative to 0 ka. The 6 ka insolation is lower over much of the Southern Hemisphere (SH) between December and June and higher between August and November. The zonal-mean difference in insolation over the whole year is plotted in Fig. 2b. There is lower insolation between 10 • N and 40 • S and higher insolation southward of 50 • S. In Fig. 2c and d, respectively, the insolation is split into two 6-month seasonal means, which coincide with the times of year when the highest insolation (October to March) and lowest insolation (April to September) occur in the SH. Between October to March (Fig. 2c), the zonal mean insolation is lower at 6 ka between 10 • N and 60 • S and higher southward of 65 • S. Conversely between April and September the insolation is higher at all latitudes between 10 • N and 90 • S (6 ka relative to 0 ka).
The only other change applied under the PMIP framework within the 6 ka simulations is a reduction in the methane concentration from 760 ppbv (0 ka) to 650 ppbv (6 ka) (Braconnot et al., 2007). Given that methane concentrations have increased from 760 ppmv (1750 CE) to approximately 1800 ppmv (April 2016) and account for approximately 17 % of the increased radiative forcing since 1750 CE (0.5 W m −2 out of approximately 3 W m −2 total; see Blasing, 2016), the climatological impact of reducing the concentration by 110 ppbv will be negligible. Finally, the calendar and seasonal definitions are the same in both the 0 and 6 ka simulations, as described in Braconnot et al. (2007).

Post-1750 CE datasets
Two other datasets are used in this study to evaluate the PMIP2 and PMIP3 experiments. The first is the Hadley Centre Sea Ice and Sea Surface Temperatures (HadISST; Rayner et al., 2003) from 1870-1899, which is instrument-based and used to evaluate the simulated sea surface temperature (SST) in the 0 and 6 ka experiments over the tropical Pacific Ocean and Southern Ocean. The period 1870-1899 is used to minimise the effect of increased SSTs from anthropogenic greenhouse gas emissions; however, data from the period 1979-2008 are also presented (see Supplement) to highlight the impact of global warming.
The second dataset is the low-level (850 hPa) zonal flow field from ERA-Interim (Dee et al., 2011) for the period 1979-2008, which is used to highlight simulated circulation errors over the tropical Pacific Ocean. As neither of the above datasets is representative of the climate at 1750 CE (as in the 0 ka simulations), they are only used to highlight known biases in the AOGCM simulations that may cause discrepancies relative to the proxy interpretations.

Analysis
To compare the model simulations with the proxy data, the area-weighted average of the climatic variable (in this case temperature or precipitation) is calculated for each simulation. An anomaly is determined by subtracting the value for 0 ka from the value for 6 ka. Two measures of multi-model agreement are then computed: i. A Student t test is used to determine whether the multi-model mean is significantly different from zero at the 5 % significance level. Each of the model simulations is assumed to be statistically independent of the others. The multi-model ensemble, however, comprises multiple different versions of the same mod-elling frameworks (examples include CCSM, CSIRO-Mk3, HadCM/GEM and MRI-CGCM; see Supplement  Table S1), and so the independence assumption may not be strictly valid. Nevertheless, each different version of the same model uses a different configuration of the parametrised physics (e.g. MRI-CGCM2.3 is configured with and without dynamic vegetation) and could be considered as a different model. The independence assumption therefore, in this situation, provides an unconditional assessment of the models' capabilities for representing the climate at 6 ka that is useful for the comparison with the available proxy data.
ii. A "model consensus" is derived by calculating the percentage of the models that agree on the sign of the temperature or precipitation anomaly (i.e. positive or negative). A value of 50 % implies that 16 models show an increase and 16 models show a decrease in temperature or precipitation at 6 ka relative to 0 ka and therefore there is no clear consensus. If the consensus is above 50 % then this indicates that ≥ 17 models agree on an increase or decrease (other examples: 21 models agree = 66 %, 25 models agree = 78 %, 29 models agree = 91 %). The consensus provides a measure of model agreement to quantify how representative the t-test result is across the model ensemble. It is also worth noting that different configurations of the same model may produce similar results and thereby increase the consensus; however, as stated above, the different physics configurations of the same model may also result in very different climatic states and so each model is treated as an independent realisation for the consensus estimate.
The two measures (above) of multi-model agreement are used to highlight whether a change in surface temperature or precipitation is a robust feature across the simulations or not. The results of the model analysis are then compared with the available proxy data. Where there is broad agreement between the models and proxies, dynamical mechanisms are presented to evaluate the proxy interpretation. Where there is disagreement between the models and proxies, possible explanations, uncertainties and biases, and the potential for future investigations are highlighted.

Temperature
Tropical north-west (TNW) and north-east (TNE) Estimates from marine sedimentary records within the Indo-Pacific Warm Pool (IPWP) suggest SSTs were broadly similar to present during 7-5 ka, with temperatures around 301-302 K (Visser et al., 2003;Stott et al., 2004Stott et al., , 2007Spooner et al., 2005;Linsley et al., 2010;Leduc et al., 2010). There is some uncertainty, however, as there is evidence of IPWP SSTs being higher (Tudhope et al., 2001 -specifically in the TNE region), equivalent (Abram et al., 2009) or lower (Gill et al., 2016) at 6 ka relative to 0 ka. In comparison, the multi-model mean differences (6 − 0 ka) in temperature for the TNW and TNE domains are −0.17 ± 0.05 and −0.21 ± 0.05 K, respectively (both statistically significant with 81 % model agreement; see Fig. 3a). The model results therefore agree with the study of Gill et al. (2016), who indicate that the lower SSTs result from an enhancement of the equatorial Pacific easterlies. Despite the apparent agreement between the PMIP models and Gill et al. (2016), there is an important caveat regarding the mean climate state over the Pacific in the AOGCMs that merits discussion. Many coupled AOGCMs are known to have poor representation of the SST field across the equatorial tropical Pacific. Typically, the SSTs are too low along the equatorial Pacific and those negative SST errors extend into the western tropical Pacific Grose et al., 2014;Irving et al., 2011;Zheng et al., 2012), which is known as the "cold-tongue" bias. Furthermore, the same errors are also visible in the PMIP2 and PMIP3 simulations used in this study (An and Choi, 2014). A simple way to remove the impact of the error is to assume that it remains unchanged in a different climatic state (such as changing the Earth's orbital parameters). Such an assumption implies that the difference between two simulated climate states is representative of the "observed" (e.g. proxy data) difference despite the initial error (as could be done when comparing to Gill et al., 2016). Nevertheless, it seems that for 6 ka conditions, the SST errors may actually be enhanced and therefore the change in the climate state is dependent on the initial error in the background state.
Previous work by Zheng et al. (2008) suggests that there was a strengthening of the Pacific trade winds in the austral spring around 6 ka, which has been attributed to a strengthening of the south-east Asian summer monsoon and is also evident in the PMIP2 and PMIP3 models An and Choi, 2014). If the tropical Pacific easterlies are already too strong in the 0 ka simulations, however, any further strengthening could enhance the existing cold-tongue bias through the Bjerknes feedback mechanism (Bjerknes, 1969;Lin, 2007). To illustrate this, the differences between the PMIP ensemble mean and HadISST (1870-1899) SSTs are plotted in Fig. 4a. Overlaid on the figure are the differences in the 850 hPa zonal wind for the ensemble mean relative to ERA-Interim (averaged over 1979. It is immediately obvious that there is a strong (< −1.5 K) SST anomaly along the western equatorial Pacific, which coincides with an easterly zonal wind bias. Recent work by Li and Xie (2014) suggests that the zonal wind error is responsible for the cold-tongue bias through the Bjerknes feedback; however, Zheng et al. (2012) suggest other mechanisms may be responsible. Regardless of the actual cause,    Figure 3. The ensemble and regional annual mean differences in (a) surface temperature (K), (b) precipitation (mm day −1 and %) and (c) 850 hPa circulation (m s −1 ) for the 6 ka simulations relative to the 0 ka simulations. In (a) blue shading (circles) indicates lower areaaveraged surface temperature and red shading (circles) indicates higher at 6 ka from PMIP (proxy) estimates. In (b), orange shading (circles) indicates lower area averaged precipitation and green shading (circles) indicates higher at 6 ka from PMIP (proxy) estimates. Grey circles indicate that proxy-derived temperature and/or precipitation at 6 ka was equivalent to 0 ka and unshaded boxes indicate changes in temperature and/or precipitation that are not statistically significant (p > 0.05) in the models. NB: in (a) the blue, grey and red boxes indicate that all three states for temperature are possible from the proxy data in the TNE at 6 ka relative to 0 ka; in (b) the orange/green rectangles denote the proxy-derived precipitation change in the northern and southern halves of the StA zone, respectively. In both (a) and (b) the values of the ensemble mean changes and the percentages of models that agree on the sign (positive or negative) of the ensemble mean temperature or precipitation differences are given. Furthermore, circles with an "X" through the middle indicate no proxy data available (both temperature and precipitation). In (c) solid and dashed contour lines indicate mean westerly and easterly flow (respectively) in the 0 ka simulations and the overlaid arrows show vector wind anomalies for 6 ka relative to 0 ka (arrow length and colour is proportional to wind anomaly strength).
both the cold-tongue bias and easterly errors are present in the PMIP ensemble in Fig. 4a. While there are uncertainties associated with both the HadISST and ERA-Interim datasets (see Rayner et al., 2003;Dee et al., 2011, respectively), the consistency between the SST (negative anomalies) and circulation (easterly anomalies) errors suggests that the biases in the models are more important than any uncertainty in the reference datasets. To this end, Fig. 4 is reproduced in the Supplement (Fig. S3a), except HadISST data are averaged over the period 1979-2008 to correspond with the ERA-Interim data. While the magnitude of the SST anomalies are larger (∼ −2.5 K), they are still evident regardless of the averaging period and the interpretation above remains valid. When the 6 ka SST and 850 hPa flow are considered relative to 0 ka, both the negative SST anomalies (relative to HadISST) and easterlies (relative to ERA-Interim) are larger in magnitude for the PMIP2/3 multi-model mean (Fig. 4b).
The difference between the 6 and 0 ka simulations (Fig. 4c) may therefore be the result of enhancing the errors that already exist in the 0 ka simulations. Hence, the models may be simulating the same conditions at 6 ka as those described in Gill et al. (2016) for the wrong reason. Conversely, if the SST at 6 ka were higher (Tudhope et al., 2001) or equivalent (Abram et al., 2009) to 0 ka then the models are not representing the IPWP correctly. It is also important to note the uncertainty in the proxy-derived SSTs in the western Pacific around 6 ka given the different estimates from Abram et al. (2009), Gill et al. (2016, and Tudhope et al. (2001). Indeed, Abram et al. (2009) show that there was a transition in the IPWP SSTs from relatively high around 6.6-6.3 ka to relatively low around 5.5 ka, compared with 0 ka. The models (given they simulate perpetual 6 ka conditions) may therefore be more representative of ∼ 5.5 ka. It should be noted that the Abram et al. (2009)    gin of the IPWP, and thus may represent a contraction of the IPWP at this specific time. Given the different proxies (e.g. corals versus sediment cores) and different localities, it is important to bring these various lines of evidence together (i.e. Abram et al., 2009;Gill et al., 2016;Tudhope et al., 2001) in order to infer the mean climate state within the TNE zone at 6 ka, such that AOGCMs can be evaluated. Overall, the discussion above indicates that fixing the cold-tongue bias is still a very high priority for the AOGCMs.
Tropical south-west (TSW) and south-east (TSE) The TSE proxy data indicate similar to present conditions at 6 ka from marine records (Bostock et al., 2006) and warmer conditions from the terrestrial island records (Woltering et al., 2014) and coral records from the Great Barrier Reef (Gagan et al., 1998(Gagan et al., , 2004, with slightly lower temperatures in the hinterlands (Haberle, 2005;Burrows et al., 2016). Speleothem records from north-west Australia (i.e. TSW; Denniston et al., 2013) suggest mean annual temperatures were equivalent to or slightly cooler at 6 ka than at 0 ka. Therefore, there appears to be a signal of lower temperatures over the land and evidence of slightly higher SSTs in the TSE at 6 ka. The models (on average) simulate lower surface temperatures at 6 ka (−0.09 ± 0.06 K) for the ensemble mean in the TSE and TSW, which are both statistically significant differences. Nevertheless, the change in temperature is very  Figure 5. The monthly, ensemble and regional mean (a) insolation (taken at 3.75 • N for insolation, black line, W m −2 ), (b) surface temperature (land and ocean combined, red line, K) and sea surface temperature (when available, amber line, K) and (c) total precipitation (blue line, mm day −1 ) and convective precipitation (turquoise line, mm day −1 ) at 0 ka (solid lines) and 6 ka (dashed lines) within the TNW box. The difference in those fields (insolation, temperature and precipitation) for 6 − 0 ka is plotted in (d).
small (< 0.1 K) and could be interpreted as similar to present (given the weak model consensus), in agreement with the proxies. It is likely that the lower 6 ka temperatures in the AOGCMs are caused by the lower annual mean insolation at these latitudes relative to 0 ka (see Fig. 2b), which agrees with the terrestrial proxies. It is clear from Fig. 4c, however, that the multi-model mean SSTs are lower in the 6 ka simulations than in the 0 ka simulations, in disagreement with the proxy evidence outlined above. Given the negative SST bi-ases in the 0 ka simulations (relative to HadISST, Fig. 4a), the model-proxy disagreement may also be related to the coldtongue bias (as described for the TNW and TNE above) and is further evidence of the need to improve the representation of the mean climate state over the Pacific.

Precipitation and circulation TNW
In the TNW, slightly drier than modern conditions are apparent at 6 ka in a lake record from Sulawesi (Russell et al., 2014); however, speleothem records from Borneo indicate similar or slightly higher annual mean precipitation at 6 ka relative to present (Partin et al., 2007;Chen et al., 2016). The multi-model mean change in precipitation for the TNW is −0.31 ± 1.34 % (Fig. 3b), which agrees with the proxy inter-pretation of similar precipitation amounts around 6 ka relative to 0 ka.
There is an interesting point raised in the work by Chen et al. (2016), that an increase in insolation over Borneo in September, October and November (SON) around 5.5 ka corresponded with increased convective activity there and also, more broadly, across the IPWP (as indicated in Partin et al., 2007). The data from the AOGCMs gives us an opportunity to test this hypothesis for 6 ka, i.e. whether the change in insolation is directly driving any changes in seasonal precipitation.
In both the 6 and 0 ka simulations, the peaks in insolation (Fig. 5a) and surface temperature (Fig. 5b) occur in September/October and March-April-May. The highest precipitation occurs in November to April (in the 0 and 6 ka simulations, Fig. 5c) for both total (blue line) and convective (turquoise lines) precipitation. The model-simulated seasonal cycle in precipitation is consistent with the observed rainfall seasonality within the TNW domain (McBride, 1998;Lee and Wang, 2014). At 6 ka (relative to 0 ka), insolation and surface temperature are higher in August to November, but precipitation is higher in November to March (Fig. 5d). Therefore, neither the mean seasonal cycle of precipitation in either period (6 and 0 ka) nor the changes in precipitation (6 ka relative to 0 ka) are driven directly by higher insolation and surface temperatures as suggested by Chen et al. (2016).
Previous work by Chang et al. (2016) and Vincent (1998) also shows that precipitation is not primarily driven by insolation in the TNW region but by the seasonal cycle in the large-scale circulation from a relatively dry southeasterly flow in April to October to a relatively moist easterly to northeasterly flow in November to March. The models also represent the seasonal change in wind direction from southeasterly during April to October (Fig. 6b) to easterlynortheasterly during November to March (Fig. 6c), which corresponds with the seasonal peak in rainfall. Furthermore, the precipitation is higher in November to March in the 6 ka simulations (relative to 0 ka), which corresponds with anomalous east to northeasterly 850 hPa flow over the TNW during October to March (Fig. 6h). It is therefore clear that the seasonality of precipitation over the TNW is driven by the large-scale circulation and not directly through higher insolation. There is one caveat, however: the change in the seasonal circulation is likely to have been caused by a strengthening of the south-east Asian monsoon in response to insolation forcing, which is also seen in the PMIP simulations (Zheng et al., 2008;An and Choi, 2014). Therefore, changes in insolation are likely to be responsible for the change in circulation plotted in Fig. 6 and thereby indirectly change the precipitation seasonality over the TNW.

TSW and TSE
Speleothem records from Flores (TSW) indicate similar amounts of precipitation fell there at 6 ka relative to 0 ka (Griffiths et al., 2009); however, marine records indicate higher precipitation in the TSW domain at 6 ka (Stott et al., 2004) as do speleothem records of warm season precipitation over north-west Australia (Denniston et al., 2013). The multi-model annual mean precipitation (Fig. 3b) is higher in the TSW at 6 ka relative to 0 ka (1.63 ± 0.92 -significant), which is primarily from an increase in October-March rainfall (5.95 ± 1.27 % -see Supplement, Fig. S1b). The models therefore largely agree with the proxy evidence.
When the TSW mean insolation and surface temperatures are plotted seasonally (Fig. 7a, b and d), the higher insolation (and surface temperatures) during October to December corresponds with higher rainfall around the same time ( Fig. 7c  and d). Furthermore, the higher simulated 6 ka rainfall is primarily from convection (turquoise line, Fig. 7c and d), indicating a thermally driven, direct response to the change in seasonal insolation at 6 ka relative to 0 ka. The models therefore agree with the proxies for higher rainfall during the warm season (i.e. October to March).
The cause of the higher rainfall at 6 ka (particularly in October to March) over TSE can be seen when the seasonal cycle of precipitation is considered (Fig. 8). The higher insolation in June to December (Fig. 8a and d) causes SST at 6 ka to be higher in August to January (Fig. 8b and d; SST response lags the insolation change), which coincides with the period where both the convective and total precipitation are higher in the 6 ka simulations relative to 0 ka ( Fig. 8c and  d). Furthermore, higher land surface temperatures also coincide with the higher convective precipitation at 6 ka relative to 0 ka (Fig. 8d), which causes an increase in low-level convergence over the land that is consistent with the stronger easterlies (see Fig. 3c and Supplement Fig. S1c). The earlier onset of the monsoon from increased continental heating, stronger onshore flow and higher SST adjacent to the land (i.e. higher evaporation) is therefore likely to be responsible for the higher precipitation over the TSE at 6 ka in the models. Conversely, the impact of reduced land and sea temperatures from April to July has little impact on precipitation during the dry season.
The precipitation characteristics for both of the TSW and TSE domains appear to respond directly to insolation (Figs. 7 and 8, respectively). In both regions, lower annual mean insolation causes surface temperatures to be lower around 6 ka relative to 0 ka; however, the lower annual mean temperatures do not result in reduced precipitation. The higher insolation from July to November at 6 ka relative to 0 ka causes the wet season precipitation to start earlier (September-October) than at 0 ka (October-November). As the response of the SST lags the insolation changes by 1-2 months, the average difference in SST in the 0 and 6 ka simulations during the middle of the wet season (December to February) is approximately ± 0.2 K (i.e. very little difference). Therefore,  Figure 7. The monthly, ensemble and regional mean (a) insolation (taken at 13.75 • S for insolation, black line, W m −2 ), (b) surface temperature (land and ocean combined, red line, K) and sea surface temperature (when available, amber line, K) and (c) total precipitation (blue line, mm day −1 ) and convective precipitation (turquoise line, mm day −1 ) at 0 ka (solid lines) and 6 ka (dashed lines) within the TSW box. The difference in those fields (insolation, temperature and precipitation) for 6 − 0 ka is plotted in (d).
given the insolation and resulting SST conditions, an overall increase in wet season precipitation occurs.

The subtropical arid zone (StA)
The StA zone incorporates much of the Australian continent and is sensitive to both the strength of the monsoon in the north and the mid-latitude westerlies in the south . Wetter conditions in the northern (tropics influenced) arid zone of Lake Eyre imply a more active monsoon at this time (6 ka relative to 0 ka), although the southern half of the arid zone is drier (Magee et al., 2004;Quigley et al., 2010;. For the whole StA region, precipitation is (on average) lower in the 6 ka PMIP  Figure 8. The monthly, ensemble and regional mean (a) insolation (taken at 13.75 • S for insolation, black line, W m −2 ), (b) surface temperature (land and ocean combined, red line, K) and sea surface temperature (when available, amber line, K) and (c) total precipitation (blue line, mm day −1 ) and convective precipitation (turquoise line, mm day −1 ) at 0 ka (solid lines) and 6 ka (dashed lines) within the TSE box. The difference in those fields (insolation, temperature and precipitation) for 6 − 0 ka is plotted in (d).
simulations by 2.20 ± 1.46 % (statistically significant) with 65 % of the models in agreement. In the StA zone north of 26.5 • S precipitation is 3.91 % lower, and only 0.4 % lower in the StA zone south of 26.5 • S. Therefore, the lower multimodel mean precipitation is primarily from the monsoondominated northern half of the StA zone, which disagrees with the relatively high precipitation from the proxies. There is, however, evidence of lower precipitation in the southern half of StA during April to September from the PMIP models in agreement with the proxies (i.e. drier winters; see . For the StA region as a whole, the April to September precipitation is 1.75 ± 1.42 % lower  Figure 9. The monthly, ensemble and regional mean insolation (taken at 23.75 • S for insolation, black line, W m −2 ), surface temperature (red line, K), total precipitation (blue line, mm day −1 ) and convective precipitation (turquoise line, mm day −1 ) at 0 ka (solid lines) and 6 ka (dashed lines) within (a) the northern half of the subtropical arid box and (b) the difference in those same fields for 6 − 0 ka for the northern half of the subtropical arid zone. Equivalent figures for the southern half of the subtropical arid zone (c and d, insolation at 28.75 • S) are also plotted. in the 6 ka simulations than in the 0 ka simulations, but the change is not significant (< 60 % of the models agree; see Fig. S2 in the Supplement). April to September precipitation is 2.35 % lower in the southern half of the StA zone (i.e. south of 26.5 • S), which is in agreement with the proxy evidence above. By looking at the seasonal characteristics of the precipitation over the northern and southern halves of the domain, as well as the arid zone as a whole, there may be evidence in the simulations to corroborate the proxy synopsis or explain the discrepancies above.
The seasonal cycles of insolation, surface temperature, convective precipitation and all precipitation are plotted for the northern arid zone in Fig. 9a. Insolation peaks in December at both 6 and 0 ka; however, surface temperatures peak in December at 6 ka and January for 0 ka. Precipitation is higher in July to December, when insolation and/or surface temperatures are higher, but there is lower rainfall from January to June, when the surface temperatures and/or insolation are lower. The precipitation in the models therefore appears to be responding primarily to the lower December to February land surface temperature (and insolation) at 6 ka (Fig. 9b). It is, however, important to consider that part of the Lake Eyre basin (which is used to infer the northern StA zone palaeoprecipitation -see above) lies within the TSE zone (see Fig. 1 in , where the multi-model mean precipitation is higher for 6 ka. It is therefore possible that the models simulate the 6 ka climate correctly, with the higher precipitation the TSE that caused the Lake Eyre basin to respond. Such a process could be evaluated in a hydrological model for the Lake Eyre basin using downscaled data from the AOGCMs but is beyond the scope of the current study and an area for future work. In the southern half of the arid zone, a similar thermally direct response to the insolation also appears in December to February, with highest monthly mean precipitation in January-February, when surface temperatures are highest (Fig. 9c). There is also a second peak in rainfall in June and July, when insolation is lowest (Fig. 9c), which is likely to be associated with extratropical systems (Catto et al., 2012). For the 6 ka multi-model mean, there is lower January to April convective precipitation at 6 ka (Fig. 9d), which is consistent with the reduced insolation and surface temperature. Conversely, the increase in insolation and surface temperature causes higher convective precipitation in October to December at 6 ka. There is also an increase in precipitation (albeit weak) in July to September, which may be indicative of an increasing influence of extratropical weather systems during the winter to early spring; however, the increase in precipitation appears to be from increased convection (turquoise line, Fig. 9d) and indicates that the higher rainfall may be a thermally direct response to the increased insolation in July to September. Overall, it appears that the lower annual mean insolation at 6 ka (relative to 0 ka) is responsible for causing lower precipitation (primarily convective) in the southern half of the StA zone and such a process (lower convective rainfall) may therefore be responsible for the drier conditions indicated by the proxies.

Temperature
Proxy records indicate that marginally higher than modern SSTs are present through the Great Australian Bight (TeS zone; Calvo et al., 2007;Lopes dos Santos et al., 2012). Terrestrial records from pollen and charcoal from the TeE and TeS, however, both indicate slightly lower temperatures at 6 ka than present (Chalson and Martin, 2012;Williams et al., 2015). The proxy-derived changes in temperature are of opposing signs over the land and ocean, however, which indicates uncertainty over whether the regional mean temperature was higher or lower at 6 ka. Lower temperature at 6 ka relative to 0 ka is simulated in the TeE (−0.11 ± 0.07 K -significant) region in agreement with the proxies. Given that the annual mean insolation is lower between 20 and 33 • S at 6 ka ( Fig. 2b) it is likely that the models and proxies are responding to that change. Conversely, in the TeS zone, the change in surface temperature at 6 ka relative to 0 ka is not significant (Fig. 3a). The annual mean difference in insolation between 6 and 0 ka is approximately zero between 33 and 40 • S, which may also explain the non-significant change in temperature from the models and uncertainty in the proxies.

Temperate east
Pollen and isotope records from North Stradbroke Island in the TeE indicate higher precipitation at 6 ka than at present , although records from Fraser Island suggest lower precipitation (Longmore, 1997;Donders et al., 2006). The pollen and charcoal records indicate drier conditions in the Sydney Basin and wetter conditions to the south at 6 ka (Chalson and Martin, 2012). The lack of a regionally coherent signal in the proxies for precipitation at 6 ka relative to 0 ka across the TeE region indicates there is uncertainty in the proxy estimates. There is also no statistically significant change in precipitation in the models (−0.46 ± 1.01 %), which is consistent with no regional consensus of higher or lower precipitation in the proxy record.

Temperate south
Sedimentology-, palaeoecology-and geochemistry-based lake records from western Victoria (TeS) indicate higher lake levels than present (Kemp et al., 2012); however, in some circumstances lower than their maximum at 7.5 ka . Furthermore, records from the western Victorian crater lakes in the TeS suggest highly variable conditions (i.e. regularly fluctuating between high and low rainfall), with a marked decrease in effective precipitation from 7 to 6 ka and around 1750 CE (e.g. Wilkins et al., 2013;Gouramanis et al., 2013). Further south, pollen and charcoal records from Tasmania reveal overall wetter conditions to the west and drier to the east (Fletcher and Thomas, 2010;Fletcher and Moreno, 2012;Jones et al., 2017). There is also evidence from both New Zealand and South America that indicates wet conditions on the western (windward) flanks of the mountains that intercept westerly flow between 40 and 44 • S across the hemisphere, suggesting that, while possibly attenuating, there was a persistence of relatively strong westerly flow at this latitude at ca. 6 ka (Fletcher and Moreno, 2012).
The strength of the westerly winds and their influence on precipitation in the TeS for both the present-day (Garreaud, 2007) and ∼ 6 ka (e.g. Fletcher and Moreno, 2012;Gouramanis et al., 2013, and references therein) periods have been widely discussed. In general, stronger (weaker) westerly flow corresponds with more (fewer) extratropical disturbances and therefore higher (lower) precipitation. Nevertheless, the correlations between precipitation and 850 hPa flow strength are ∼ 0.3 over southern Australia (Garreaud, 2007), which indicates that other processes must be important. Indeed, Garreaud (2007) show that the correlations are higher when both the 850 hPa zonal flow and relative humidity are considered over Australia, which indicates that moisture availability is also an important limiting factor for precipitation over the TeS and not just the circulation strength. Both Fletcher and Moreno (2012) and Gouramanis et al. (2013) indicate that precipitation was higher around 6 ka but for different reasons. Fletcher and Moreno (2012) state that it is "a phase of enhanced westerly flow that led to a moisture increase in western Tasmania and south-west Victoria", whereas Gouramanis et al. (2013) state that "strengthening of this narrow, warm surface water [Leeuwin Current] may increase convection and precipitation locally, resulting in shifted pressure gradients causing wetter conditions in southeast Australia. . . and south-west Australia" (i.e. not directly caused by stronger westerlies). The AOGCMs on the other hand indicate no change/slightly higher precipitation (only +0.3 %) at 6 ka relative to 0 ka with 65 % of models hav- ing higher precipitation, despite there being weaker westerly flow (see Fig. 3c). A similar anti-correlation exists in October to March where 71 % of the models simulate higher precipitation at 6 ka (+2.14 ± 1.39 % -significant, Supplement Fig. S1b) even though the westerlies are weaker (Supplement Fig. S1c). Interestingly, there is little consensus on reduced April to September rainfall (56 % of the models simulate reduced precipitation, Supplement Fig. S2b) when the westerlies should have their greatest influence on southern Australia precipitation through frontal cyclones (Catto et al., 2012). As the physical processes responsible for precipitation can be diagnosed directly from the PMIP simulations, there is an opportunity to explain (i) why the mid-latitude westerlies are likely to have been weaker at 6 ka relative to 0 ka and, (ii) despite those weaker westerlies, precipitation over the TeS domain may have been equivalent to or slightly higher at 6 ka than around 0 ka -in agreement with both Fletcher and Moreno (2012) and Gouramanis et al. (2013).
Explaining the above may provide a useful alternative mechanism to account for periods in the past when the rainfall-westerly wind strength relationship may weaken or break down.
In order to explain the cause of the weaker mid-latitude westerlies at 6 ka (relative to 0 ka), the thermal wind balance equation is considered: where U g / z is the change in the zonal geostrophic wind ( U g , m s −1 ) with height ( z, m -also called the vertical shear of the geostrophic wind with respect to height), f is the Coriolis parameter (s −1 ), T is the temperature at a reference point (K), g is the acceleration due to gravity (m s −2 ), and T / y is the change in surface temperature (K) per distance of latitude (m). Following Eq. (1), reducing equatorial surface temperatures and increasing them at high latitudes would reduce the T / y term (i.e. weaker Equator-to-pole temperature gradient). If all other parameters are held fixed then U g / z will also reduce. In order to identify whether the Equator-to-pole temperature gradient has changed, the following calculation is undertaken: where T 30N to 30S is the area averaged surface temperature (K) between 30 • N and 30 • S, T 60S to 90S is the area averaged surface temperature (K) between 60 and 90 • S and DT ep is the Equator-to-pole temperature difference (K). The values of DT ep are calculated for each individual model and plotted in Fig. 10 for the 0 ka simulations (white box, left axis), the 6 ka simulations (amber box, left axis) and the difference in DT ep for 6 ka relative and 0 ka (pink box, right axis). The median DT ep from the models is 44.7 K for 0 ka and 44.0 K for 6 ka; however, all 32 models simulate a reduction in the difference in temperature between the Equator and poles (pink box plot -upper whisker is less than zero). The weaker Equator-pole temperature gradient is consistent with lower insolation in the tropics and higher insolation at high latitudes at 6 ka compared to 0 ka (Fig. 2b). Therefore, given the insolation and surface temperature characteristics of the PMIP simulations visible in Figs. 2 and 10 (respectively), the  Figure 11. The monthly, ensemble and regional mean (a) insolation (taken at 36.25 • S for insolation, black line, W m −2 ), (b) surface temperature (land and ocean combined, red line, K) and sea surface temperature (when available, amber line, K) and, (c) total precipitation (blue line, mm day −1 ) and convective precipitation (turquoise line, mm day −1 ) at 0 ka (solid lines) and 6 ka (dashed lines) within the TeS box. The difference in those fields (insolation, temperature and precipitation) for 6 − 0 ka is plotted in (d).
westerly winds should be weaker in the mid-latitudes at 6 ka relative to 0 ka (as seen in Fig. 3c).
So, given that there is a physically plausible reason why the westerlies would have been weaker at 6 ka relative to 0 ka, why is there little change to the annual mean rainfall (or even a tendency for a small increase)? To answer this question, the seasonal cycles of insolation, temperature and precipitation are plotted for the 0 and 6 ka simulations in Fig. 11a-c and for 6 − 0 ka in Fig. 11d. At 0 and 6 ka the insolation peaks in December and is lowest in June; however, at 6 ka insolation is higher in July to November and lower in December to May than at 0 ka. There is also a shift in the seasonal surface temperature and local SST with higher surface temperatures at 6 ka relative to 0 ka between August and January. Between August and November there is lower convective precipitation and, between December and June, higher convective precipitation in the 6 ka simulations relative to 0 ka. This suggests that the non-convective rainfall (e.g. frontal rain) is not changing. Furthermore, in April to June, the convective precipitation has increased more than the total precipitation change, which indicates that the non-convective rainfall has actually reduced during those months. So, there has been an increase in convective rainfall in the PMIP models (through increased frequency and/or intensity), which has compensated for any reduction in precipitation caused by the weaker westerlies and is consistent with Gouramanis et al. (2013). Nevertheless, the flow in the models is still predominantly westerly over the TeS at 6 ka (in agreement with Fletcher and Moreno, 2012) despite being slightly weaker than at 0 ka. Overall, the model assessment given above provides an example of where the relationship between westerly wind strength and precipitation could weaken through changes in convective precipitation (which can be directly diagnosed from the models). Such a process would be very difficult to decipher from the proxies and shows where model data may be useful to elucidate the key processes that induce past climatic states.
There is, however, an important caveat associated with the model-derived precipitation estimates for 6 ka. The coarse resolution of the models means that surface topographical features on the land are not represented well. Therefore, the impact of such topography on the prevailing circulation and precipitation would also be misrepresented. Areas where such a problem may be important are over the Great Dividing Range and Tasmania. It is logical to conclude that the misrepresentation of topography may also be contributing to any interpretation of the models' simulated climate. The only way to resolve such an issue would be to run high-resolution regional climate model simulations over the TeS zone driven by output from fully coupled, multi-millennial transient global model simulations. A measure of the time-dependent change in the circulation and its interaction with the land surface could then be assessed. Such model simulations have been shown to improve the representation of present day precipitation over Southern Alps of New Zealand (Ackerley et al., 2012) and have also been applied to simulations of 6 ka .

North Southern Ocean (NSO)
The proxies suggest that SSTs were around 284.2 K in the NSO region (Barrows et al., 2007) for the annual mean at 6 ka. The HadISST-derived 1870-1899 SSTs are 284.0 K for the NSO domain. Therefore, SSTs in the NSO are almost identical at 6 and 0 ka. The models also simulate little change in SST between 6 and 0 ka (0.04 ± 0.04; see Fig. 3a). Therefore, the models and proxies agree that the SSTs in the NSO region at 6 ka are likely to have been very similar to 0 ka. The lack of any SST change in NSO is also consistent with the insolation changes between 40 and 50 • S, which are negligible (see Fig. 2b).

South Southern Ocean (SSO)
SSTs in the SSO zone for February are estimated to have been approximately 278.7 K (Crosta et al., 2004). The HadISST-derived 1870-1899 mean February SSTs are also 278.7 K for the SSO region. Therefore, as with the NSO, SSTs in the SSO are almost identical at 6 and 0 ka from the proxy evidence.
The February SSO multi-model mean SSTs in the 0 and 6 ka simulations are 280.3 and 280.4 K, which are > 1.5 K higher than the HadISST and Crosta et al. (2004) estimates given above. Furthermore, the models also simulate higher SSTs in February at 6 ka relative to 0 ka (approximately 0.11 K) with 63 % agreement. Higher SSTs are also visible for the annual mean (0.20 ± 0.07 K -significant; see Fig. 3a), with very high model agreement (91 %). It appears that the model SSTs are responding to the higher annual mean insolation at 6 ka relative to 0 ka (see Fig. 2b), which is not seen in the proxy estimate. While there are acknowledged spatial and temporal gaps in the proxy data for the Southern Ocean that may cause the slight disagreement outlined above , there are also significant known deficiencies in the model simulations within this region.
Trenberth and Fasullo (2010) have shown that the cloud cover fraction over the Southern Ocean is too low within the CMIP3 models, which leads to a positive bias in the amount of solar radiation absorbed at the ocean surface. Furthermore, this cloud-related bias in the absorbed solar radiation (approximately +8 W m −2 ; Schneider and Reusch, 2016;Flato et al., 2013) is still present in the next generation of coupled climate models (CMIP5 and PMIP3) and corresponds with positive SST biases (see Fig. 9.2 and 9.5 in Flato et al., 2013). Approximately 91 % (70 %) of the models simulate higher February (annual) mean SSTs between 50 and 60 • S at 0 ka relative to HadISST, which is consistent with the known cloud cover and radiation errors described above. Other work by Wang et al. (2014) shows that the positive SST biases are strongest in the SH summer and autumn, which corresponds with the time that the cloud cover related errors in the absorbed solar radiation are at a maximum (i.e. when insolation is highest). Nevertheless, Wang et al. (2014) attribute the Southern Ocean warm bias to errors in the meridional overturning circulation and not the cloud radiative errors.
Overall, it is clear that there are large errors in the Southern Ocean circulation within AOGCMs that could be caused by different processes (e.g. cloud radiative forcing and the meridional overturning circulation), which may enhance (or even dampen) the SST response to a change in insolation. Therefore, while the simulated higher 6 ka SST (relative to 0 ka) in the SSO for February (∼ 0.1 K) and annually (∼ 0.2 K) may be a manifestation of these errors (especially given the large mean state bias). Given limited proxy records and the fact that they are only representative of summer conditions, estimates of seasonal or annual mean SSTs at 6 ka (and other periods) are necessary to enable better model-proxy validation. Nevertheless, improving the representation of the atmosphere and ocean at high southern latitudes should be a priority given the known errors that exist in both CMIP3 and CMIP5 (Flato et al., 2013;Trenberth and Fasullo, 2010).

Future directions
Given the assessment above, this section focusses on some key opportunities for future work from both the proxy and modelling communities in order to provide a better platform to undertake fully integrated studies.

Proxies
Due to the sampling resolution of most of the proxy records and the response time of the systems from which they come, it is very difficult to reconcile seasonal variability. Exceptions to this are tree rings, coral and speleothem records, although their coverage within the vast Australasian region is sparse. One area for improvement in proxy reconstructions is a clear understanding of the season that is represented by (particularly) the biological archives, e.g. the season of pollen production and dispersal or invertebrate blooms, as is already being undertaken as part of the PAGES 2k initiative (PAGES 2k Consortium, 2017). In many cases this may be known for the organisms in question, but often not adequately described in the reconstructions, or the ranges not considered. There is great potential to re-interrogate the proxy records in view of the model outputs with regard to changes in the seasonality of the signal. There are also possibilities of deriving more quantitative data from proxies, either through the use of transfer functions or calibrating geochemical variability on the organisms directly, to look more at seasonal variability. Suitable proxies for these studies include tree rings, speleothems and molluscs that show clear incremental growth in addition to coral records. As always, more robust chronologies can only benefit high-resolution palaeoclimatic work in addition to targeting areas of climatic sensitivity and poor geographic coverage.

Models
The original work by Reeves et al. (2013a) compares the regional surface temperature and effective precipitation characteristics of one period relative to a previous one and not relative to the present day. Therefore, the first logical step would be to run time slice simulations of each of the time periods discussed in Reeves et al. (2013a). A more ambitious plan would be to develop transient model simulations of the last 35 kyr, which would allow a direct comparison with the OZ-INTIMATE synthesis; however, although feasible, such simulations would be computationally very expensive. Despite that, there are some multi-millennial model simulations that have already been undertaken (e.g. Wagner et al., 2007;Menviel et al., 2011;Liu et al., 2009) and have also been planned as part of PMIP4 (Ivanovic et al., 2016), which may provide a template for other modelling groups to follow. Such simulations would also be useful to investigate decadal to millennial variation in large-scale climate modes (such as El Niño-Southern Oscillation).
In order to acquire high-resolution model data to compare with the proxies, higher-resolution global AOGCMs (e.g. Fallah et al., 2016) or regional climate models (RCMs, dynamical downscaling) could be employed (e.g. Ackerley et al., 2013;Wagner et al., 2012). This may be particularly important over complex terrain (such as Tasmania and the Great Dividing Range) where this study has identified differences between the proxy reconstructions and modelled climates. Nonetheless, bias corrections need to be applied to the boundary data from the global AOGCM used to drive the RCM (not done in Ackerley et al., 2013;Fallah et al., 2016;Wagner et al., 2012), otherwise the simulations may only reproduce the existing systematic model errors, but at higher resolution.
Future model developments should also aim to include proxy system models (PSMs; Dee et al., 2016) or use available model variables to produce pseudoproxies Smerdon, 2012). Both approaches provide a more integrated way of comparing AOGCMs and proxies rather than the relatively simple comparative study employed here. Nevertheless, the usefulness of both PSMs and pesudoproxies will also depend on any systematic biases inherent to the AOGCM being used, as well as on any uncertainties inherent in the formulation of the PSMs and pseudoproxies themselves. Therefore studies that validate climate models against proxies directly (as done here) are still a necessity.

Summary and conclusions
This study aimed to investigate the AOGCM-simulated (from the PMIP ensemble) climate state within the geographical regions defined by Reeves et al. (2013a) relative to the available proxy data for the mid-Holocene (6 ka) within those regions. Where the models and proxies agreed, the influence of the external forcing (insolation) or circulation (atmospheric dynamics) was presented in order to evaluate the proxy interpretation. Where there was uncertainty associated with the model simulations and/or proxies (e.g. the cold-tongue bias in the models and a lack of consensus in the proxy estimates), the reasons for this were discussed in order to highlight opportunities for further research.
The main results of this study are as follows: -In most of these areas, surface temperature and precipitation respond directly to the changes in insolation. The one exception was the tropical north-west (TNW), where precipitation was driven by circulation change and not directly from the insolation.
-The simulated change in climate at 6 ka is sensitive to the "cold-tongue bias" in the tropical Pacific appar-ent in the 0 ka simulations. It is the enhanced easterly flow over the tropical Pacific (from the stronger southeast Asian monsoon) that enhances the error. However, complexities in the comparison of the multiple proxy records also need to be considered.
-Annual mean rainfall over the temperate south (TeS) appears to be unchanged for 6 ka relative to 0 ka despite weaker westerly flow. Higher convective precipitation balances a reduction in precipitation from extratropical systems. When modern-day analogues of relating precipitation to westerly wind strength break down, the models may offer a useful alternative mechanism (e.g. changes to convective precipitation). Conversely, the coarse resolution of the AOGCMs may mean they are not representing the climate in topographically diverse regions (e.g. Tasmania and the Great Dividing Range).
-Southern Ocean SSTs are higher at 6 ka relative to 0 ka from the increased insolation; however, there is no evidence from the proxy data for this. The discrepancy may be due to the poor model representation of clouds over the Southern Ocean and/or the proxy reconstruction only being representative of February conditions (i.e. not the annual mean).
Overall, this study shows that "upscaling" proxy reconstructions (to provide coherent regional information) is the most direct method to compare with coarse-resolution climate models. Where there is agreement, the models can be used to verify the inferred dynamical processes from the proxies -especially for determining the role of external forcing mechanisms. Where there is disagreement, the proxies and models can be evaluated further to understand (and hopefully address) the cause of the mismatch, which may relate to inherent uncertainties within the proxies of regionally idiosyncratic responses to incident atmospheric circulation (see Goodwin et al., 2013;Lorrey et al., 2008Lorrey et al., , 2014. It is also important to compare regionally coherent proxy information against similar regions within an AOGCM. Comparing single proxy reconstructions with individual model grid points is unlikely to yield useful results as the models cannot simulate the climate at such scales (grid spacings of > 100 km). Finally, there is also a clear need for acquiring seasonally resolved proxies to evaluate the impacts of orbital variations on seasonality, which the models are capable of simulating.